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ABSTRACT 


Although  the  progress  made  in  the  past  thirty  years 
to  describe  the  transient  nature  of  gas  flow  through  porous 
media  has  been  somewhat  staggering,  very  little  experimental 
data  exist  to  substantiate  the  validity  of  the  results  ob¬ 
tained.  Tests  were  performed  on  a  number  of  different  cores 
in  the  presence  of  a  residual  fluid  saturation  and  the  results 
analyzed  on  the  basis  of  an  effective  porosity  and  permeabil¬ 
ity  to  gas.  Originally  it  was  intended  that  the  presence  of 
the  immobile  fluid  would  simulate  a  connate  water  saturation. 
Kerosene,  nonetheless,  was  chosen  for  this  second  phase  in 

order  to  prevent  damaging  the  core  by  the  expansion  of  in- 

(20  33 ) 

herent  clay  materials  on  contact  with  water  '  .  The 

theory,  however,  shows  that  the  method  of  analysis  for  a 
residual  oil  phase  is  also  applicable  to  the  presence  of  a 
residual  water  saturation. 

The  validity  of  the  dimensionless  cumulative  pro¬ 
duction  function  and  the  radius  of  drainage  equation  derived 
for  the  constant  terminal  pressure  case  was  tested  by  com¬ 
parison  with  experimental  data  obtained  for  heterogeneous 
porous  systems.  It  was  found  that  the  presence  of  hetero¬ 
geneities  seriously  limited  the  applicability  of  the  common 
dimensionless  parameters  used  in  normal  reservoir  calculations. 
A  numerical  solution  was  also  effected  and  the  results  compared 


to  experimental  data.  The  large  discrepancies  obtained  were 
attributed  mainly  to  limitations  in  the  experimental  appara¬ 
tus.  The  presence  of  the  fluid  saturation  further  accentu¬ 
ated  these  limitations  by  increasing  the  reaction  time  of 
the  equipment  to  changes  in  pressure. 
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INTRODUCTION 


Because  of  the  complexity  of  fluid  flow  through  con 
solidated  porous  media,  original  investigations  into  the 
transient  phenomena  yielded  results  purely  of  an  academic  in¬ 
terest.  It  was  not  until  1949,  when  Van  Everdingen  and 
(43 ) 

Hurst  published  their  classical  solutions  for  slightly 

compressible  fluids,  that  a  practical  method  was  available 
to  analyze  common  reservoir  engineering  problems.  Their  pre¬ 
sentation  of  the  dimensionless  cumulative  production  and  pres 
sure  functions  greatly  facilitated  the  prediction  of  reser¬ 
voir  histories.  Later  other  investigators  extended  this  work 
to  include  the  flow  of  gases.  However,  even  before  these 
solutions  are  applied  to  a  practical  problem,  further  simpli¬ 
fying  assumptions  must  be  made  with  regard  to  the  physical 
properties  of  the  porous  medium.  Heterogeneities  are  usually 
combined  into  single  effective  values  of  porosity  and  per¬ 
meability.  For  gas  reservoirs  connate  water  saturations  are 
considered  only  to  reduce  the  effective  porosity  available 
to  the  gaseous  phase,  and  the  effects  of  solubility  are  com¬ 
pletely  neglected. 

The  advent  of  high  speed  electronic  computers  facil 
itated  a  numerical  approach  for  a  direct  solution  to  the  non¬ 
linear  gas  flow  equation.  The  introduction  of  finite  differ¬ 
ence  approximations  for  the  derivatives  enabled  the  analysis 
of  more  sophisticated  problems.  It  now  became  possible  to 


■ 
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incorporate  such  pressure  dependent  terms  as  gas  viscosity 
and  gas  compressibility  factors  into  the  final  solution. 

The  phenomenon  of  slip  and  its  influence  on  the  transient 
response  was  analyzed.  Heterogeneities  in  the  form  of  per¬ 
meability  distributions  were  introduced  and  the  range  of 
analyses  extended  to  two  and  three  dimensions. 

Although  these  concepts  have  been  the  basis  for 
numerous  engineering  studies,  to  date  very  little  experi¬ 
mental  data  have  been  provided  to  substantiate  the  applica¬ 
tion  of  the  developed  theory  to  practical  situations.  Thus, 
the  prime  objective  of  this  investigation  was  to  obtain 
transient  data  on  a  number  of  different  reservoir  cores,  and 
to  determine  how  well  the  theory  would  describe  the  experi¬ 
mental  results  obtained.  A  second  immobile  phase  was  intro¬ 
duced  to  simulate  a  connate  water  saturation,  and  attempts 
were  made  to  analyze  the  theory  proposed  to  describe  the  con¬ 
stant  terminal  pressure  case.  In  particular,  emphasis  was 
placed  on  the  dimensionless  cumulative  production  function 
and  the  radius  of  drainage  equation  derived  for  an  infinite 
linear  system. 
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LITERATURE  REVIEW 

A  review  of  the  literature  was  conducted  in  order 
to  determine  the  extent  of  the  theoretical  and  experimental 
techniques  available  to  describe  the  transient  phenomenon. 

The  survey  revealed  that  the  entire  analysis  of  transient 
isothermal  flow  of  gas  through  porous  media  is  based  on  the 
solution  of  the  general  diffusivity  equation  and  its  inherent 
assumptions 

k  9  P 

V  ( —  PVP)  =  4  —  (-)  (1) 

yZ  9 1  Z 

where 

k  =  effective  permeability  of  the  porous  media, 
darcies 

y  =  viscosity  of  the  gas,  centipoises 
Z  =  compressibility  factor  of  the  gas 

4  =  effective  porosity  of  the  porous  media, 

fraction 

P  =  pressure,  atmospheres 

t  =  time,  seconds 

with  all  lengths  measured  in  centimeters. 

A  historical  review  of  the  solution  to  this  problem 
may  be  divided  into  two  parts,  the  point  of  separation 
roughly  occurring  with  the  advent  of  high  speed  electronic 
computers . 
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Early  Empirical  and  Analytical  Investigations 

Assuming  that  the  properties  of  the  porous  media 
and  flowing  gas  are  constant,-  equation  (1)  may  be  simplified 
to 

2  2 

v  P Z  =  c  —  (2) 

8 1 

where  C  is  a  constant  dependent  upon  the  properties  of  the 
gas  and  the  porous  media0  Because  of  the  non-linear  charac¬ 
teristics  of  equations  (1)  and  ( 2 )  ,  analytical  solutions  to 

these  relationships  have  not  yet  been  obtained. 

( 34 )  • 

Muskat  and  Botset  in  1931  made  one  of  the  first 

attempts  at  solving  this  transient  behavior  of  gases  in  porous 
media.  By  way  of  a  dimensional  analysis  of  the  problem,  they 
obtained  the  general  law  for  the  flow  of  gas  through  a  sand 
as 

6P2  =  k(pv)n  (3) 

where  pv  is  the  mass  velocity,  p  being  the  density,  and  v 

the  effective  velocity0  The  left  hand  side  of  the  equation 

is  equivalent  to  the  difference  in  pressures  squared 
2  2 

(P^  -  P 2  )  and  n  is  an  empirical  constant  determining  the 
nature  of  flow,  varying  between  the  limits  1  and  2,  and 
corresponding  to  completely  viscous  and  completely  turbu¬ 
lent  flow  respectively.  The  final  solution  for  a  radial 
system  was  obtained  by  applying  a  continuous  succession  of 


- 
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steady  states,  each  governed  by  the  equation 

=  0  (4) 

where  S  represents  the  space  variable,  the  time  being  in¬ 
troduced  as  a  parameter  through  the  boundary  conditions. 

This  approximation  is  equivalent  to  assuming  that  the  velo¬ 
city  of  flow  is  negligible  when  compared  to  the  propagation 

of  disturbances  through  the  porous  material.  In  a  later 

( 35 ) 

analysis  Muskat  introduced  the  concept  of  a  "radius  of 

drainage",  which  is  considered  to  recede  to  the  boundary 
of  the  reservoir  as  the  steady  state  condition  is  estab¬ 
lished.  In  this  instance  no  gas  is  removed  from  any  point 
within  the  system  until  the  radius  of  drainage  has  passed 

through  that  point.  Similar  attempts  were  made  to  apply 

(36) 

these  concepts  to  two  phase  flowv  with  a  minor  degree  of 
success . 

( 19 ) 

Hetherington  et  al  developed  empirical  equa- 

.  i 

tions  relating  cumulative  production  and  future  production 
rates  from  gas  reservoirs  as  a  function  of  time  for  a  linear 
system  with  a  sealed  external  boundary.  Again  the  approxi¬ 
mation  of  a  continuous  succession  of  steady  states  was  em¬ 
ployed,  with  further  assumptions  being  that  the  average  pres¬ 
sure  over  the  entire  system  is  negligibly  different  from  the 
pressure  at  the  external  boundary,  and  the  downstream  pres¬ 
sure  is  small.  However,  unlike  the  work  of  their  predeces- 


div 


3P 

( - ) 

3  S 


1/n 
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sors,  the  equations  advanced  by  Hetherington  et  al  are  based 
mainly  on  the  characteristic  back-pressure  equation. 


Q 

-  =  a  (P 


Pf2)N 


(5) 


Q  =  total  volume  of  gas  produced  at  time  t 
Pe  =  pressure  at  the  external  boundary  of  the 
system 


f 

a 


=  pressure  at  the  outlet  side  of  the  system 
=  experimental  constant  describing  the 
geometry  of  the  system 


N  =  experimental  exponent  describing  the 
nature  of  the  flow 

and  as  a  consequence  are  void  of  any  complicated  integration 
or  series  terms.  The  validity  of  these  equations  was  favor¬ 
ably  tested  by  comparison  with  experimental  results  and  the 
concept  of  a  "steady  state  decline"  introduced.  That  is  to 
say,  upon  the  inception  of  production  there  would  initially 
exist  a  period  of  transient  conditions  referred  to  as  "the 
recession  of  the  radius  of  drainage",  followed  by  a  longer 
period  of  steady  decline  in  which  the  rate  of  change  of 
pressure  with  respect  to  time  would  be  independent  of  posi¬ 
tion  in  the  system. 

The  classical  solutions  of  Van  Everdingen  and 

(8,25,43)  r  ,  .  ,  r\ 

Hurst  for  equation  (6), 


9  P  1  9P 
9  r  2 


4>  w  c  9P 


r  9r 


k  9 1 


(6) 
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which  describes  the  transient  radial  flow  of  a  slightly 
compressible  fluid  through  porous  media,  have  been  the 
foundation  for  the  analysis  of  many  reservoir  engineering 
problems  concerning  water  influx  and  fluid  flow  (both  oil 
and  gas)  into  a  well  bore*  A  slightly  compressible  fluid 
has  been  defined  as  that  which  is  governed  by  the  equation 

where  c  is  the  compressibility  of  the  fluid  and  the  sub¬ 
script,  o,  refers  to  some  standard  condition. 

Upon  further  examination  it  is  noted  that  equation 
(6)  is  linear  and  as  such  is  amenable  to  analytic  treatment. 

Although  explicit  solutions  for  equation  (6)  had  been  ob- 

(22  35) 

tained  in  prior  investigations  v  '  ,  these  were  much  too 

involved  and  cumbersome  to  be  of  any  value  in  handling  prac¬ 
tical  engineering  problems.  Basically  the  approach  of  Van 
Everdingen  and  Hurst  consisted  of  converting  equation  (6) 
and  the  prescribed  boundary  conditions  into  an  ordinary  dif¬ 
ferential  equation  by  applying  the  Laplace  transformation. 
The  transformed  equation  was  then  solved  analytically  and 
the  result  converted  to  "real"  time  by  way  of  Mellin's  in¬ 
version  formula.  This  analysis  involved  the  introduction 
of  two  dimensionless  parameters  as  functions  of  dimension¬ 
less  time:  Q  ,  the  dimensionless  cumulative  influx  term; 

tD 

and  P  ,  the  dimensionless  pressure  term,  where 
tD 


1  +  c (P  - 


Po> 


(7) 


-  8 


t_  i— 


'D 


r 

J 

o 


Q, 


9P 


^  9  (~L 

rw 


dt 


D 


r=r 


'D 


(Pf  -  Pe> 


(8) 


and 


r  =  radial  distance 


w 


'D 


=  inner  radius 


=  dimensionless  time 


=  constant  flowing  pressure  at  the  inner  radius 
=  the  initial  reservoir  pressure 


p.  =  f(t  ) 
D 


(9) 


the  exact  nature  of  the  function  on  the  right  hand  side  of 
equation  (9)  being  dependent  upon  the  value  of  t  itself 
and  the  boundary  conditions  applied.  The  final  results  were 
presented  as  working  tables  and  charts  of  dimensionless  para¬ 
meters,  thus  eliminating  the  time  consuming  process  of  evalu¬ 
ating  complicated  series  solutions  involving  Bessel  functions 
It  was  soon  discovered  that  the  Van  Everdingen  and 
Hurst  solutions  could  be  applied  to  gas  flow  through  porous 
media  if  equation  (2)  is  linearized  in  pressure  squared  by 

introducing  an  average  constant  reservoir  pressure  for  the 

(12  13  25) 

time  interval  considered'  '  f  .  Rewriting  equation  (2) 


provides 


2  2 

VP  =  C’ 


9P 


9 1 


(10) 


* 
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where  C'  incorporates  the  assumed  average  pressure,  P.  The 

2 

net  effect  is  to  replace  the  pressure,  P ,  by  P  and  the  fluid 
compressibility,  c,  by  the  inverse  of  P  in  all  the  relation¬ 
ships  developed  for  a  slightly  compressible  fluid.  The 
method  of  solution  now  lends  itself  to  a  trial  and  error  pro¬ 
cedure  in  an  attempt  to  choose  an  appropriate  value  for  P. 

Finite  Difference  Approximations  and  Numerical  Solutions 

(1) 

Aronofsky  and  Jenkins  obtained  numerical  solu¬ 
tions  to  equation  (2)  for  linear  systems  of  infinite  and 
finite  length,  in  which  the  initial  and  terminal  pressures 
and/or  rates  were  specified.  The  non-linear  gas  flow  equa¬ 
tion  was  expressed  in  finite  difference  form,  and  the  pres¬ 
sure  at  time  t  +  At  determined  explicitly.  In  order  to  in¬ 
sure  numerical  stability  it  was  necessary  to  limit  the  size 
of  the  time  increment  as  determined  by  equation  (11). 

a  2 

At  <.  -  ( Ax ;  (11) 

4 

2  <f>y 

a  =  -  ,  where  P.  is  the  initial  pressure  in  the  system. 

P  .k *  1 

i 

A  comparison  of  the  results  obtained  with  an  analytical  re¬ 
lationship  for  transient  fluid  flow  revealed  the  necessity 
for  specifying  the  pressures  involved  in  order  to  uniquely 
define  the  solution  for  the  flow  of  gases.  The  basis  of  this 
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work  was  later  extended  to  include  transient  radial  flow 

(2  23 ) 

through  porous  media  '  .  However,  in  this  instance  the 

restrictions  imposed  on  At  for  a  predictive  solution  were 

partially  overcome  by  the  introduction  of  a  fictitious  well 

bore  radius.  Steady  state  conditions  were  assumed  to  exist 

within  the  bounds  of  this  radius  and  it  was  necessary  to 

correct  the  final  numerical  solutions  to  account  for  the  gas 

produced  from  this  region. 

The  work  of  Roberts ^ 40 ^  is  significant  in  that  he 
2 

showed  that  9P  /9x  is  a  measure  of  the  flux,  and  the  non¬ 
linear  effect  of  the  gas  flow  equation  acts  in  such  a  way 
as  to  retard  the  rate  at  which  the  pressure  drops  within 
the  system.  A  solution  was  accomplished  by  applying  a  step 

by  step  linearization  process  to  the  diffusivity  equation 

(18) 

and  integrating  the  result  numerically.  Green  and  Wilts v 
solved  the  transient  problem  for  gas  flow  by  satisfying  the 
finite  difference  equation  with  an  electrical  analogy.  The 
limitations  of  this  approach  are  obvious. 

The  classical  numerical  analysis  of  equation  (2) 
for  both  the  linear  and  radial  systems  was  performed  by 
Bruce,  Peaceman,  Rachford,  and  Rice ^  in  1953.  An  implicit 
procedure  was  proposed  which  greatly  relaxes  the  stringent 
control  on  the  maximum  size  of  At.  This,  however,  necessi¬ 
tates  the  solution  of  a  tri-diagonal  matrix  with  each  itera¬ 
tion  of  the  numerical  procedure. 
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Cornell  and  Katz  developed  a  graphical  method 

of  solution  in  an  endeavor  to  account  for  turbulence ^ ^ ^  in 
the  vicinity  of  the  well  bore.  It  was  assumed  that  steady 
state  turbulent  flow  exists  in  the  region  near  the  well  bore 
and  that  laminar  transient  flow  prevails  at  some  larger  dis¬ 
tance  away  from  the  well.  The  procedure  consists  of  first 
solving  for  the  unsteady  state  case  and  then  correcting  the 

result  for  turbulence. 

( 3 ) 

Aronofsky  extended  the  scope  of  his  earlier 

numerical  approach  to  include  the  effects  of  molecular 

streaming  or  slip.  It  had  been  shown  in  a  prior  investiga- 

(27) 

tion  by  Klinkenberg  that  the  effective  permeability  to 

gas  is  a  linear  function  of  the  mean  pressure  in  the  system 
and  is  dependent  on  the  nature  of  the  flowing  gas. 


b 

k  =  K  ( 1  +  — ) 


(12) 


b  = 


P  = 


effective  permeability  to  gas 
absolute  (liquid)  permeability  and  is 
dependent  on  the  geometry  of  the  porous 
media  only 

a  constant  dependent  on  the  nature  of  the 
gas  and  the  porous  material 
mean  pressure 


The  modus  operandi  consisted  of  incorporating  equation  (12) 
in  the  finite  difference  approximations  and  solving  the 
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resultant  expression.  It  was  shown  from  the  results,  and 
later  by  experimentation'^5^,  that  the  slip  effect  is  negli¬ 
gible  for  pressures  normally  encountered  in  gas  reservoirs 

even  though  deviations  in  the  laboratory  may  be  significant. 
(28) 

Kolada  J  combined  the  effects  of  molecular  streaming  and 

turbulence  in  his  steady  state  analysis  of  the  viscous  and 

( 9 ) 

visco-inertial  regions  of  flow.  Chwyl  then  extended  this 
concept  to  the  transient  case  by  combining  equation  (12)  with 
the  well-known  Forchheimer  equation.  His  numerical  solution 
showed  that  the  effects  of  turbulence  tend  to  counteract 
the  effects  of  slip. 

A  similar  technique  to  account  for  pressure  depen¬ 
dent  variables  was  employed  by  Aronofsky  and  co-authors  to 

( 4 ) 

determine  the  non-ideal  behavior  of  real  gases  in  linear 
( 5 ) 

and  radial  systems.  Both  the  compressibility  factor  Z 
and  the  viscosity  y  were  assumed  to  be  simple  functions  of 
pressure,  and  in  particular 


1  -  mP 

(13) 

d  +  nP 

(14) 

where  m,  d,  and  n  are  constants  dependent  on  the  nature  of 
the  gas. 

(16) 

Eilerts  in  1964  combined  the  major  contribu¬ 

tions  of  previous  investigators , in  particular  Bruce  et  al, 
and  developed  a  numerical  method  for  the  solution  of 
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equation  (1)  for  a  linear  system.  Effect  was  given  to  the 
pressure  dependency  of  the  compressibility  factor,  the  vis¬ 
cosity,  and  the  effective  gas  permeability,  and  to  the  change 
of  absolute  permeability  with  position.  The  pressure  de¬ 
pendent  properties  were  combined  into  one  simple  function  of 

pressure  before  expressing  the  differential  equation  in  finite 

( 7 ) 

difference  form.  Carter  extended  this  analysis  to  two- 
dimensional  steady  state  flow  and  effected  a  solution  by  the 
method  of  Saul'ev,  commonly  referred  to  as  ADEP  (Alternating 
Direction  Explicit  Procedure).  Like  Eilerts,  Carter  combined 
the  pressure  dependent  properties  into  one  term,  which  he  de¬ 
fined  as 

P  AdX 

*(P)  =  /  -  (15) 

Pc  u(A)ZU) 

P  in  equation  (15)  is  any  fixed  value  of  pressure. 

Although  the  above  review  is  by  no  means  complete, 
it,  nonetheless,  serves  to  indicate  the  degree  of  sophistica¬ 
tion  to  which  the  transient  analysis  of  gas  flow  has  pro¬ 
gressed.  Of  major  concern,  however,  is  the  apparent  lack 
of  experimental  evidence  to  support  many  of  the  concepts 
proposed . 
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THEORY 


Basic  Diffusivity  Equation  for  Gas  Flow 

Before  an  analysis  of  the  basic  diffusivity  equa¬ 
tion  for  gas  flow  can  be  effected,  it  is  imperative  to 
examine  the  limitations  imposed  on  this  relationship  by  way 
of  its  derivation.  A  review  of  the  basic  laws  describing 
the  flow  of  fluids  reveals  that  equation  (1)  is  derived  by 
combining  the  equation  of  continuity  governing  the  flow  of 
fluids  through  porous  media  of  constant  porosity, 

-  3P 

V  •  (pv)  =  -<J>  -  ( 16 ) 

at 

the  equation  of  state  for  a  gas, 

MP 

P  =  -  (17) 

ZRT 

and  Darcy's  Law  neglecting  gravitational  effects. 

k 

v  =  -  -  VP  (18) 

u 

R,  in  equation  (17),  is  the  universal  gas  constant.  M,  the 
molecular  weight  of  the  gas,  and  T,  the  temperature,  are 
assumed  to  be  constant. 

To  be  strictly  correct  the  empirical  relationship 
expressed  by  equation  (18)  should  be  replaced  by  the  classical 
hydrodynamic  equations  of  Navier-Stokes .  However,  the  solu¬ 
tion  of  this  system  of  equations  for  a  geometry  as  complex 
as  a  porous  medium  would  prove  to  be  insurmountable. 
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If  the  permeability,  the  viscosity,  and  the  com¬ 
pressibility  factor  are  assumed  to  be  constant,  equation  (1) 
can  be  re-written  as 

<j>y  3P 

V  •  (PVP)  = -  (19) 

k  at 

Applying  the  mathematical  relationship  expressed  by 

VP2  =  2PVP  (20) 

equation  (19)  can  be  simplified  to 

9  9  2  4>y  3P 

VZPZ  = -  (21) 

k  at 

Finally,  multiplying  the  right  hand  side  of  equation  (21) 

by  P/P  and  assuming  an  average  pressure, P  in  the  denomina- 

2 

tor,  yields  an  equation  linear  in  P  . 

9  9  4>y  3P2 

V  P Z  =  —  -  (22) 

kp  at 

i 

It  can  now  be  shown  that  the  study  of  transient 
gas  flow  in  the  presence  of  a  residual  fluid  saturation  is 
theoretically  substantiated  by  the  expressions  describing 
multiphase  flow  (Appendix  I).  The  method  of  obtaining  a 
solution  is  the  same  as  for  a  single  phase  system.  However, 
now  the  compressibility  of  the  single  phase  is  replaced  by 
the  total  compressibility  of  the  system,  and  the  mobility 
is  replaced  by  the  total  mobility  of  the  flowing  fluids. 

For  the  case  of  gas  flowing  in  the  presence  of  a  residual 
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fluid  saturation,  only  the  mobility  and  compressibility  of 
the  gaseous  phase  need  be  considered.  The  latter  approxima¬ 
tion  is  attributed  to  the  small  values  of  liquid  compres¬ 
sibilities  when  compared  to  the  compressibilities  of  gases. 


Dimensionless  Cumulative  Production 


Following  the  procedure  outlined  by  Van  Everdingen 
(43) 

and  Hurst  for  the  constant  terminal  pressure  case,  it 

can  be  shown  that  for  linear  gas  flow 


n 


c 


cj>LA 


P 


ZfRT 


-  P  )Q 
e 


(23) 


where 

fined 


is 


the  dimensionless  cumulative  production  de- 

8P 

tn  *3 (x/L) 'x=0 

Qt  =  JD  -  dt  (24) 

D  o  (Pf  -  Pe) 

=  cumulative  production  in  moles 
=  length  in  inches 

=  average  pressure  from  the  dimensionless 
time  expression  in  psia 
=  area  in  inches 

=  constant  flowing  pressure  at  the  downstream 
end  in  psia 

=  constant  pressure  at  the  external  boundary 
in  psia 
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R  =  universal  gas  constant 
T  =  temperature  in  °R 
<p  =  porosity,  fraction 
Z ^  =  compressibility  factor  at 


Radius  of  Drainage 

(44) 

Van  Poolen's  investigation  of  the  various 

equations  describing  the  radius  of  drainage  and  stabiliza¬ 
tion  time  reveals  that  the  final  expressions  obtained  are 
dependent  upon  the  criteria  used  to  define  both  the  radius 
of  drainage  and  stabilization  time.  Consider  the  solution 
for  the  constant  terminal  pressure  case  governing  the  linear 
flow  of  gas  in  an  infinite  porous  system. 


where 


2  2 
P  (x,t)  -  P0 


erf  c 


(25) 


kPt 


u<px 


(26) 


Defining  the  radius  of  drainage  as  the  point  where 


P2 (x,t)  -  Pe2 


=  0.01 


f  e 

and  following  a  procedure  similar  to  Jones 
equation  for  the  radius  of  drainage  becomes 


(24) 


(27) 


,  the  final 


4  /ktPyp  <j> 


(28) 
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and  upon  re-arranging  and  solving  for  t,  the  stabilization 
time  becomes 


16kP 

The  units  of  equations  (28)  and  (29)  are  in  the  c-g-s  system. 
It  should  be  noted  that  the  derivation  of  equation  (28)  in¬ 
volves  making  the  assumption 


1.82  -  2.0  (30) 

Consequently,  a  more  accurate  solution  to  the  radius  of  drain¬ 
age  would  be 

xd  =  3.6  4/ktF/uicf)  (31) 


Numerical  Solution 

Assuming  a  constant  compressibility  factor  and 
defining  the  compressibility  of  the  gas  as  1/P,  equation  (1) 
for  a  linear  system  becomes 

a  k  (x)  9P2  <p  ap2 

—  - )  =  =  -  (32) 

9x  V (P)  9x  P  at 

where  the  permeability  is  now  a  function  of  position  and  the 
viscosity  a  function  of  pressure.  A  finite  difference  approxi¬ 
mation  was  obtained  for  equation  (32)  by  a  method  similar  to 
(16) 

Eilerts v  7  (Appendix  II),  and  a  computer  program  written  to 
numerically  solve  the  result  (Appendix  III)  for  a  number  of 
boundary  conditions. 
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EXPERIMENTAL  APPARATUS 

Description  of  Cores 

Tests  were  performed  on  a  total  of  four  different 
cores,  an  Indiana  limestone,  two  Berea  sandstones , and  a 
"Series"  core  composed  of  a  Berea  sandstone  and  a  section  of 
Alundum.  Excluding  the  Series  core  for  the  moment,  the  re¬ 
maining  cores  were  set  in  5-inch  I.D.  oil  well  production 
casing,  using  an  epoxy  resin  to  fill  the  annular  space. 

The  method  for  pouring  the  resin  is  described  in  detail  by 
Mackett'^^.  Once  a  permanent  set  had  been  achieved,  five 
1/4-inch  diameter  pressure  taps  were  machined  at  convenient 
intervals  along  the  length  of  the  core.  The  depths  of  these 
pressure  taps  extended  through  the  casing  and  epoxy  resin, 
and  approximately  1/4  inch  into  the  core  material  itself. 

The  ends  of  the  core  samples  were  fractured,  as  opposed  to 
machining,  to  prevent  damage  to  the  core  faces. 

The  core  holder  consisted  of  two  heavy  duty  mild 
steel  flanges,  specially  designed  to  accommodate  a  metal 
insert  at  the  inner  faces.  These  flanges  were  held  together 
on  either  side  of  the  core  by  eight  3/4-inch  diameter  mild 
steel  rods  threaded  at  both  ends.  For  instances  where  the 
core  face  protruded  past  the  length  of  the  casing,  spacer 
rings  were  inserted  between  the  flange  and  the  casing.  Metal 
to  metal  contact  was  avoided  by  using  Garlock  rubber  gaskets, 
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which  provided  an  excellent  seal  to  pressures  as  high  as 
900  psig.  Centrally  located  pressure  taps  permitted  the 
flow  of  gas  through  the  flanges,  and  a  second  pressure  tap 
provided  a  direct  means  of  measuring  the  pressure  at  the 
core  faces. 

A  Series  core,  shown  in  Figure  1,  was  assembled 
to  study  the  effects  of  heterogeneities  on  the  transient 
behavior  of  a  porous  system.  The  Berea  sandstone  was  intro¬ 
duced  in  two  distinct  sections  in  an  attempt  to  determine 
the  influence,  if  any,  of  the  Berea-Alundum  interface.  The 
three  sections  of  core  were  assembled  to  simulate  two  dif¬ 
ferent  permeabilities  in  series  and  set  in  4-inch  I.D.  pro¬ 
duction  casing.  Five  3/8-inch  diameter  pressure  taps  were 
then  drilled  along  the  length  of  the  core  to  facilitate  the 
measuring  of  pressure  distributions.  Unlike  the  previous 
arrangements,  two  inches  of  thread  were  machined  at  both 
ends  of  the  casing  to  accommodate  metal  caps.  These  pro¬ 
vided  the  final  seal  for  the  core  faces.  It  should  be  noted 
that  each  of  these  metal  caps  contained  only  a  single  cen¬ 
trally  located  pressure  tap  to  allow  for  communication  be¬ 
tween  the  core  and  the  remaining  portion  of  the  experimental 
apparatus.  Consequently  it  was  necessary  to  install  a  "T" 
at  both  ends  of  the  Series  core  arrangement  to  provide  a 
means  for  obtaining  the  absolute  pressures  at  these  points. 
Two  complete  sets  of  experimental  data  were  obtained  for 
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Figure  1  Cross  Section  of  Series  Core  Assembly 


22 


this  core;  one  with  the  Alundum  section  upstream  of  the 
Berea  sandstone  (AUS) ,  and  a  second  with  the  core  positions 
reversed  (ADS) . 

Measuring  Devices 

The  transient  behavior  of  commercially  pure  nitro¬ 
gen  gas  was  measured  on  a  model  906C  Honeywell  "Visicorder" , 
which  provided  a  permanent  record  of  the  pressure  drop  across 
each  section  of  the  core  as  a  function  of  time  (refer  to 
Figure  2) .  Five  ±15  psid  Statham  PM80TC  differential  pres¬ 
sure  transducers  provided  the  electrical  input  into  the 
Visicorder.  A  sixth  ±50  psid  Statham  PM60TC  transducer  was 
required  at  the  producing  end  of  the  core,  because  of  the 
large  pressure  drops  experienced  in  this  region.  A  0-18v 
Kepco  Laboratories  power  supply  provided  the  electrical  in¬ 
put  into  a  central  panel,  which  in  turn  controlled  the  in¬ 
dividual  inputs  into  each  transducer.  Inputs  into  the 
transducers  were  monitored  on  a  Non-Linear  Systems  Model 
481  digital  voltmeter. 

The  temperature  of  the  flowing  gas  stream  was 
sensed  by  two  iron-constantan  thermocouples  and  the  electri¬ 
cal  output  measured  on  a  Leeds  and  Northrup  millivolt 
potentiometer.  As  shown  in  Figure  2,  the  upstream  and  down¬ 
stream  pressures  at  the  core  faces  were  measured  with  bour¬ 
don  tube  pressure  gauges.  A  0-1500  psig  pressure  regulator 
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Schematic  Diagram  of  Experimental  Apparatus 


inserted  at  the  nitrogen  supply  permitted  direct  control  of 
the  flowing  upstream  pressure.  More  precise  control  was 
effected  by  means  of  a  continuously  bleeding  0-50  psig 
Moore  "Nullmatic"  pressure  regulator  installed  downstream 
of  the  surge  tank.  All  pressure  lines  consisted  of  1/8-inch 
O.D.  Autoclave  stainless  steel  tubing  and  fittings,  rated  at 
10,000  psig.  Finally,  a  Toggle  valve  was  introduced  on  the 
producing  side  of  the  core.  The  snap-open  characteristics 
of  this  valve  enabled  the  simulation  of  the  boundary  condi¬ 
tion  prescribing  the  sudden  decrease  in  pressure  at  the 
downstream  location. 
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EXPERIMENTAL  PROCEDURE 


Effective  Porosity 


The  effective  porosity  of  each  core  sample  was  de¬ 


termined  directly  by  expanding  a  gas  isothermally  into  a 
known  volume.  This  necessitated  the  construction  of  a 
special  porosimeter  which  is  shown  in  Figure  3.  The  extran¬ 
eous  pressure  taps  in  each  core  were  sealed  by  utilizing  a 
combination  of  steel  ball  bearings  and  female  fittings.  The 
procedure  consisted  of  closing  valve  A  and  introducing  a 
controlled  amount  of  nitrogen  into  the  system  by  means  of  a 
pressure  regulator.  Several  minutes  were  usually  required 
to  attain  stabilized  conditions,  and  the  equilibrium  pres¬ 
sure  was  read  from  a  bourdon  tube  pressure  gauge.  A  vacuum 
was  simultaneously  drawn  on  the  expansion  chamber  and  the 
result  measured  on  a  mercury  manometer.  Valve  A  was  then 
opened  allowing  direct  communication  between  the  core  sample 
and  expansion  chamber.  Once  equilibrium  conditions  had  been 
achieved,  the  final  pressure  was  measured  on  the  manometer 
and  the  porosity  calculated  from 


where 


porosity  in  percent 


if  , 

■ 


0-200  psig 
Bourdon  Tube  Gauge 
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Figure  3  Schematic  Diagram  of  Porosimeter 
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V 


SL 


N 


=  bulk  volume  of  the  core  obtained  from 
caliper  measurements 

V  =  volume  of  the  expansion  chamber  and  inter¬ 
connecting  pressure  lines,  including  the 
manometer 

=  increase  in  volume  of  the  total  system  due 
to  the  change  in  the  level  of  mercury  in  the 
manometer 

=  additional  apparent  void  volume  of  the  core 
due  to  the  nitrogen  source  lines,  core 
holder  spacer  rings,  etc. 

=  initial  pressure  in  the  core  sample 
=  initial  pressure  in  the  expansion  chamber 
=  final  equilibrium  pressure  in  the  system 
Equation  (33),  as  presented,  is  valid  for  any  system  of  con¬ 
sistent  units. 

Several  runs  were  performed  on  each  core  at  various 
pressure  levels  and  the  results  from  equation  (33)  averaged 
to  obtain  a  final  porosity.  All  barometric  pressure  readings 
were  corrected  for  latitude  and  temperature  throughout  the 
entire  course  of  the  investigation.  It  was  noted  that  by 
following  the  above  procedure,  initial  nitrogen  pressures  of 
less  than  170  psig  provided  final  equilibrium  pressures  below 
atmospheric.  Consequently  no  vacuums  were  drawn  on  the  ex¬ 
pansion  chamber  for  these  runs  in  order  to  prevent  damage 
from  occurring  to  the  pressure  gauge.  Equation  (33),  however, 
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is  still  valid. 

Saturating  the  Core 

In  order  to  obtain  the  liquid  or  absolute  perme¬ 
ability  of  each  core,  it  was  first  necessary  to  saturate  the 
entire  porous  network  with  kerosene .  This  was  accomplished 
by  completely  evacuating  the  core  with  a  vacuum  pump;  then 
opening  a  valve  to  allow  the  kerosene  to  naturally  imbibe 
into  the  pore  spaces,  all  the  while  continuously  drawing  a 
vacuum  at  the  opposite  end  of  the  sample.  The  extent  of  the 
vacuum  drawn  was  measured  on  a  mercury  manometer  as  shown 
in  Figure  2.  The  above  procedure  was  continued  until  break¬ 
through  of  kerosene  was  detected  in  the  vacuum  flask G  Ap¬ 
proximately  40-50  psig  pressure  was  externally  applied  to 
the  surface  of  the  kerosene  in  the  lucite  container  to  aid 
in  the  imbibition  process  of  the  tighter  cores.  A  simple 
before-and-after  weighing  of  the  lucite  container  and  vacuum 
flask  permitted  a  material  balance  to  be  performed.  This, 
along  with  a  knowledge  of  the  density  of  kerosene,  provided 
a  means  for  calculating  the  volume  of  fluid  imbibed.  Next 
the  line  was  broken  downstream  of  the  lucite  cylinder  and  a 
WII  Ruska  Proportionating  Pump  introduced.  The  core  was  then 
flushed  with  approximately  450  cc  of  kerosene,  the  efflux 
discharged  into  the  vacuum  flask,  and  a  second  material  ba¬ 
lance  performed  to  account  for  any  volumetric  discrepancies. 


■V 
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Finally  kerosene  was  "squeezed"  into  the  core  at  150  psig. 

The  difference  between  the  initial  and  final  pump  readings 
at  0  psig  represented  a  further  fluid  loss  to  the  core.  By 
combining  the  results  of  the  above  procedure  and  accounting 
for  the  presence  of  kerosene  in  the  flow  lines,  etc.,  it  was 
possible  to  obtain  a  check  on  the  effective  porosity  previously 
determined.  In  all  cases  agreement  was  obtained  to  within 
1.2  percent  of  the  total  porosity. 

Absolute  Permeability 

To  determine  the  absolute  permeability,  kerosene 
was  injected  into  the  respective  cores  at  some  known  preset 
rate  utilizing  the  Ruska  pump.  The  upstream  and  downstream 
temperatures  and  pressures  at  the  steady  state  were  recorded 
by  means  of  the  thermocouples  and  pressure  gauges  shown  in 
Figure  2,  and  the  permeability  calculated  from  the  inte¬ 
grated  form  of  Darcy's  Law 

qyL 

K  =  -  (34) 

AAP 

where  AP  is  the  pressure  drop  equal  to  (P2  -  P^) .  A  periodic 
check  of  was  sufficient  to  determine  as  to  when  steady 
state  conditions  had  been  attained.  Several  runs  were  per¬ 
formed  at  different  injection  rates,  and  the  results  of 
equation  (34)  averaged  to  obtain  a  final  value  for  the  ab¬ 
solute  permeability.  An  experimentally  determined  correla- 
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tion  provided  a  means  for  obtaining  the  viscosity  of  kerosene 
at  the  measured  temperatures.  The  length  and  cross-sectional 
area  of  each  core  were  determined  by  direct  measurement.  In 
cases  where  V  ^  exceeded  250  psig,  the  upstream  gauge  was  re¬ 
placed  by  a  0-1500  psig  Heise  pressure  gauge. 


Relative  Permeability  Ratio  and  Residual  Fluid  Saturation 

The  composite  relative  permeability  ratio  of  each 
core  was  determined  by  the  familiar  "modified  external  gas 
drive  technique",  and  the  results  interpreted  according  to 
the  procedure  outlined  by  Welge^^,  and  based  on  his 
classical  relationship 


(f  •  ,  ) 

oil 


(35) 


where 


S 

av 


S 


2 


(f 


oil 


A  typical  set  of 


average  gas  saturation  in  the  core 
gas  saturation  at  the  producing  end  of 
the  core 

fraction  of  oil  flowing  at  the  producing 
end  of  the  core 

cumulative  gas  injected  in  pore  volumes 
at  the  mean  pressure  in  the  core 
results  are  shown  in  Figure  32,  Appendix 


IV.  Nitrogen  gas  was  used  to  displace  the  kerosene  from 
the  core,  and  a  gravitational  separation  of  the  two  phases 
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effected  by  means  of  the  apparatus  shown  in  Figure  4.  The 
oil  production  was  determined  directly  from  the  scale  read¬ 
ings  on  the  inverted  burettes  and  the  cumulative  gas  injected 
measured  on  a  calibrated  Sargent  Wet  Test  Meter.  The  pro¬ 
cess  was  continued  until  the  liquid  phase  completely  ceased 
to  flow.  Because  of  the  large  volumes  of  nitrogen  injected, 
it  was  necessary  to  saturate  the  gas  by  bubbling  through 
water  in  order  to  preserve  the  calibration  of  the  wet  test  meter. 
Finally  a  volumetric  balance  on  the  kerosene,  accounting  for 
liquid  within  the  lines,  etc.,  provided  a  means  for  deter¬ 
mining  the  residual  fluid  saturation.  In  all  cases  attempts 
were  made  to  use  pressures  higher  than  those  expected  in 
ensuing  tests,  in  order  to  eliminate  the  possibility  of  addi¬ 
tional  fluid  being  produced. 

Effective  Gas  Permeability 

The  effective  gas  permeability  of  each  core  in 
the  presence  of  the  residual  fluid  saturation  was  determined 
from  a  standard  Klinkenberg  Plot,  utilizing  the  integrated 
form  of  Darcy's  Law  for  the  flow  of  gases. 


T  Z 
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(36) 


T 
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T 


temperature  of  the  gas  in  the  core 
temperature  at  which  gas  is  measured 
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av 
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(p2-Pi> 
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A 


=  pressure  at  which  gas  is  measured 
=  average  pressure  in  the  core  equal  to 
(P2  +  P^/2 

=  average  gas  compressibility  factor  based 

P  and  T 
av 

=  gas  compressibility  factor  at  P^  and  T 

=  pressure  drop  across  the  core 

=  measured  gas  flow  rate  at  T  and  P 

3  mm 

=  gas  viscosity  at  P  and  T 
^  1  av 

=  cross-sectional  area  of  the  core 


L  =  length  of  the  core 

However,  for  equation  (36)  to  be  valid,  the  flow  must  be  vis¬ 
cous.  Consequently  efforts  were  made  to  determine  whether 
or  not  the  flow  rate  data  obtained  fell  in  the  viscous 
region.  The  regions  of  viscous  and  visco-mertial  flow  were 


established  by  the  method  of  Dranchuk  and  Kolada 


(15) 


The 


standard  back-pressure  plot  in  Figure  5  indicates  that  the 
region  of  flow  substantiating  the  use  of  equation  (36)  is 
very  small.  A  re-evaluation  of  the  theory,  however,  suggests 
that  a  modified  back-pressure  plot  accounting  for  the  depen¬ 
dence  of  Z,  y,  and  k  on  pressure  would  provide  a  more  accurate 
criterion  for  establishing  the  turbulent  region.  The  method 
now  lends  itself  to  a  trial  and  error  procedure  in  an  attempt 
to  establish  the  values  of  b  and  K  in  equation  (12) .  Usually 
the  lower  flow  rates  of  Figure  5  are  used  to  provide  an  initial 
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Figure  5 


Back-Pressure  Plot  -  Cut  Berea  Sandstone 
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guess  for  b  and  K  (see  Figure  6).  The  final  results  of  the 
above  procedure  are  shown  in  Figures  6  and  7.  Note  that  in 
Figure  7  only  the  highest  two  flow  rates  deviate  from  the 
45°  line  defining  the  viscous  region.  Consequently  equation 
(36)  is  valid  for  all  but  these  last  two  flow  rates.  The 
effective  permeability  to  gas  can  then  be  obtained  by  deter¬ 
mining  the  value  for  the  inverse  of  the  mean  pressure  and 
reading  the  result  directly  from  Figure  6.  The  above  analy¬ 
sis  was  performed  on  the  Indiana  limestone,  the  "Cut"* Berea 
sandstone,  and  the  high  permeability  Berea  sandstone.  Gas 
permeabilities  for  the  remaining  cores  were  determined  from 
equation  (36)  for  single  point  flow  tests. 

Stabilization  pressures  and  temperatures  at  the 
steady  state  were  determined  for  eight  different  flow  rates 
using  the  apparatus  shown  in  Figure  2.  The  rate  of  gas 
flow  was  measured  using  a  triple  element  Model  20ER"Vol-o- 
Flow"  meter  calibrated  for  the  range  0-290  cfh. 

Calibration  of  Transducers 

The  method  of  calibrating  a  ±15  psid  transducer 
consisted  of  exerting  a  14.0  psi  pressure  differential  across 

*  The  significance  of  the  adjective  Cut  is  explained  in 
the  section  entitled  "Analysis  of  Data  and  Discussion 
of  Results" . 
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Figure  6  Klinkenberg  Plot  -  Cut  Berea  Sandstone 
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Figure  7  Modified  Back-Pressure  Plot  -  Cut  Berea  Sandstone 
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the  element,  and  adjusting  the  input  voltage  to  give  a  full 
scale  deflection  on  the  Visicorder.  The  input  was  monitored 
on  the  digital  voltmeter  and  the  result  recorded  for  future 
reference.  The  pressure  differential  was  then  decremented 
by  1  psi  at  a  time  and  the  corresponding  deflections  noted. 

A  plot  of  the  deflection  versus  the  pressure  differential 
provided  the  final  calibration.  The  calibration  of  the  ±50 
psid  transducer  varied  only  in  that  a  50.0  psi  pressure  dif¬ 
ferential  was  used  to  obtain  a  full  scale  deflection,  and 
subsequent  readings  taken  at  pressure  decrements  of  5  psi. 
All  transducers  were  calibrated  against  a  0-150  psig  Seegers 
Standard  Precision  Test  Gauge.  The  element  of  the  gauge  was 
precise  to  within  0.25  psig. 

Testing  Procedure 

Since  the  apparatus  shown  in  Figure  2  was  located 
in  a  constant  temperature  cabinet,  the  commencement  of  each 
day's  testing  was  preceded  by  a  30-minute  warm  up  period, 
required  for  the  electrical  apparatus  and  experimental  en¬ 
vironment  to  stabilize.  The  transducers  were  then  isolated 
and  the  inputs  individually  reset  to  give  a  full  scale  de¬ 
flection  for  the  previously  mentioned  pressure  differentials. 

For  each  experimental  run  the  initial  condition 


P (x,o) 


P . 
l 


(37) 
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was  simulated  by  pressuring  the  core  with  nitrogen  to  some 
controlled  level.  Equilibrium  conditions  were  indicated  by 
an  equalization  of  upstream  and  downstream  pressures,  and  by 
noting  a  zero  deflection  of  all  galvonometers  on  the  Visi- 
corder.  Next  the  Visicorder  itself  was  adjusted  for  the 
desired  paper  speed  and  the  vertical  line  timer  set.  One  or 
two  trial  runs  were  usually  sufficient  to  determine  which 
settings  produced  the  best  results.  The  run  was  then  initi¬ 
ated  by  simultaneously  tripping  open  the  Toggle  valve  and 
setting  the  paper  control  mechanism  into  motion.  During  the 
course  of  each  run  it  was  usually  necessary  to  change  the 
paper  speed  or  the  vertical  line  timer  set  or  both  at  least 
once.  Such  changes  were  all  duly  noted  to  facilitate  future 
analysis . 

The  condition  specifying  a  constant  pressure  at 
the  external  boundary  was  achieved  by  means  of  the  Nullmatic 
regulator.  The  sealed  boundary  was  obtained  by  merely  clos¬ 
ing  the  valve  downstream  of  this  regulator  prior  to  initiat¬ 
ing  a  run.  In  order  to  effect  an  analysis  of  the  data,  the 
actual  value  of  the  pressure  was  required  at  one  point  or 
more  in  the  core  for  all  values  of  time.  The  constant  up¬ 
stream  pressure  was  sufficient  in  the  former  case.  However, 
for  the  sealed  boundary  the  variation  of  either  the  upstream 
or  downstream  pressures  was  required  as  a  function  of  time. 
Both  were  obtained  by  periodic  observations  of  the  two 
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pressure  gauges  and  the  actual  time  recorded  from  a  stop 
watch  activated  at  the  beginning  of  each  run.  It  was  not 
possible  to  maintain  a  constant  pressure  at  the  producing 
boundary  for  the  high  permeability  cores.  Consequently,  it 
was  necessary  to  repeat  the  same  run  at  least  once,  in  order 
to  accurately  establish  the  variation  in  upstream  and  down¬ 
stream  pressures  with  time  for  the  case  of  the  sealed  bound¬ 
ary.  Although  only  one  of  these  pressure  histories  was  re¬ 
quired,  the  second  served  as  a  check  for  the  first.  A  simi¬ 
lar  history  of  the  variation  in  downstream  pressure  with 
time  provided  a  check  for  the  constant  pressure  case.  The 
actual  constant  terminal  pressure  case  was  obtained  only  for 
the  tight  cores.  In  such  instances  the  pressure  at  the  pro¬ 
ducing  end  stabilized  at  0  psig  and  remained  there  for  the 
duration  of  the  run. 

Because  of  the  limitations  imposed  by  the  trans¬ 
ducers,  pressures  in  the  range  of  50  psig  were  used  through¬ 
out  the  course  of  this  investigation.  The  flowing  gas  tem¬ 
perature  was  measured  with  the  thermocouples  shown  in  Figure 
2,  and  the  cabinet  temperature  obtained  from  a  mercury  ther¬ 


mometer  . 
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ANALYSIS  OF  DATA  AND  DISCUSSION  OF  RESULTS 

Experimental  Data 

As  previously  mentioned,  the  Visicorder  charts 
provided  a  permanent  record  of  the  galvonometer  deflections 
as  a  function  of  time  on  light  sensitive  paper.  These  de¬ 
flections  corresponded  to  the  pressure  differentials  sensed  by 
the  transducers  for  the  respective  sections  of  the  core.  To 
determine  values  of  pressure  it  was  first  necessary  to  con¬ 
vert  the  deflections  to  pressure  drops  by  means  of  the  cali¬ 
bration  charts  obtained  earlier.  Thus,  knowing  the  pressure 
drop  from  section  to  section  within  the  core,  the  pressure 
distribution  was  obtained  directly  by  the  summation  of  these 
differentials  starting  at  a  point  of  known  pressure.  The 
above  procedure  was  repeated  for  a  number  of  different 
times  and  a  graphical  representation  of  the  transient  be¬ 
havior  obtained  for  each  core.  A  typical  set  of  results 
for  the  cases  of  a  constant  pressure  at  the  external  boundary 
and  a  sealed  external  boundary  are  shown  in  Figures  8  and  9 
respectively.  Similar  results  for  the  remaining  cores  are 
presented  in  Appendix  IV.  Deflections  were  read  to  within 
0.1  of  a  division.  This  provided  a  resolution  of  0.1  psid 
for  the  PM80TC  transducers  and  a  resolution  of  approximately 
0.4  psid  on  the  PM60TC  transducer  .  However,  because  of  the 
procedure  of  obtaining  pressure,  an  error  in  measurement  for 
one  transducer  was  passed  on  to  subsequent  pressure  values. 


Constant  Pressure  at  the  External  Boundary 
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Figure  8  Transient  Behavior  -  Berea  Sandstone 


Time  t  (min.)  Sealed  External  Boundary  p.  =  50.3  psig 
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Figure  9  Transient  Behavior  -  Berea  Sandstone 
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For  example,  if  the  PM80TC  transducers  were  all  in  error  by 
0.1  psid  and  the  PM60TC  transducer  was  in  error  by  0.4  psid, 
the  final  pressure  reading  determined  could  be  in  error  by  as 
much  as  0.9  psi. 

A  further  inspection  of  Figures  8  and  9  indicates 
a  drastic  deviation  from  the  transient  results  expected  for 
a  homogeneous  system.  Originally  this  was  attributed  to  the 
failure  of  totally  eliminating  the  mechanical  end  effect  due 
to  machining  of  the  core  faces.  As  a  result,  the  core  casing 
and  epoxy  were  cut  and  the  core  itself  sheared  at  a  distance 
of  one  pressure  tap  in  from  each  end,  in  an  attempt  to  once 
again  eliminate  these  discontinuities.  Subsequent  testing 
of  this  "Cut"  Berea  sandstone  indicated  that  although  the 
overall  effective  permeability  had  increased,  discontinuities 
were  still  prevalent  at  the  ends  (see  Figures  29a  and  29b, 
Appendix  IV).  Thus,  it  was  concluded  that  the  heterogeneous 
nature  of  the  core  was  either  a  natural  phenomenon  or  a  di¬ 
rect  result  of  damage  from  prior  testing. 

In  all  instances  tests  with  a  constant  pressure 
at  the  external  boundary  were  run  until  steady  state  condi¬ 
tions  were  achieved.  The  steady  state  pressure  profile  then 
provided  a  means  for  determining  the  effective  permeability 
of  each  section  of  the  respective  cores.  Where  possible, 
flow  rates  at  the  steady  state  were  obtained  by  direct  mea¬ 
surement  using  the  Vol-o-flow  meter.  However,  in  a  number 
of  instances  an  adjustment  to  the  equipment  prevented 
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measurements  of  the  flow  rates  from  being  obtained.  For  such 


cases  the  upstream  and  downstream  pressures  were  used  to 

2  2 

evaluate  (P^  -  )  and  the  value  of  the  desired  flow  rate 

read  from  a  back-pressure  plot  similar  to  Figure  5.  Thus, 
knowing  the  flow  rate  through  each  section  of  the  core  and 
the  corresponding  pressure  drop,  the  permeability  distribu¬ 
tion  was  determined  by  utilizing  equation  (36) .  The  results 
for  each  core  are  shown  in  Table  1.  The  last  column  of  Table 
1  represents  the  effective  gas  permeability  of  core  as  de¬ 
termined  from  equation  (38),  which  is  derived  for  permeabili¬ 
ties  in  series. 


Lef  f  ective 


n 

I  L. 

i=l  - 
n  L 

l  t 

0  T  J\ 
1=1  : 


(38) 


n  =  number  of  sections 
L  =  length  of  each  section 
k  =  permeability  of  each  section 
These  overall  core  permeabilities,  shown  in  the  last  column 
of  Table  1,  were  then  averaged  with  the  gas  permeabilities 
determined  as  described  in  the  section  titled  "Experimental 
Procedure".  This  provided  a  suitable  value  for  the  charac¬ 
teristic  gas  permeability  of  each  core.  Other  properties 
and  characteristics  of  the  cores  tested  are  listed  in 


Table  2. 
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Table  1 


Permeability  Variations  in  Test  Cores 


Overall  Core 


Core 

*Section 

(number) 

Length 
(in. ) 

Permeability 

(md) 

Permeability 

(md) 

Berea 

1 

4^ 

1.73 

Sandstone  1 

2 

6 

98.3 

3 

6 

97.3 

4 

6 

22.3 

5 

6 

8.21 

6 

1.49 

4.59 

Cut  Berea 

1 

3-7/16 

2.52 

Sandstone 

2 

3^ 

15.4 

3 

3% 

38.8 

4 

3h 

37.8 

5 

3h 

86.7 

6 

3h 

7.21 

9.11 

Berea 

1 

3h 

387 

Sandstone  2 

2 

3h 

484 

3 

3h 

526 

4 

3h 

536 

5 

3h 

431 

6 

3-15/16 

515 

474 

Series  Core 

1 

8-5/8 

315 

AUS 

2 

9 

458 

3 

9 

663 

4 

9 

838 

5 

9 

1039 

6 

9-3/8 

962 

602 

Indiana 

1 

3-11/16 

1.76 

Limestone 

2 

3h 

5.20 

3 

3h 

127 

4 

3h 

253 

# 

5 

3h 

1.93 

6 

3-10/16 

2.72 

3.58 

*  Section  numbers  begin 

from  the 

downstream  end 

of  the  core. 
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Experimental  Results 

Initially  the  tests  were  designed  to  study  the 
applicability  of  the  theory  developed  for  transient  gas 
flow  through  a  homogeneous  linear  porous  system  with  a 
liquid  saturation.  However,  all  of  the  cores  tested  proved 
to  be  of  a  heterogeneous  nature.  In  fact,  the  Series  core, 
the  only  core  assembled  with  a  deliberate  discontinuity , 
provided  results  more  characteristic  of  a  homogeneous  system 
when  the  Alundum  section  was  located  in  the  upstream  posi¬ 
tion,  Figures  24  and  25  of  Appendix  IV  for  the  Series  core 
AUS  show  that  there  are  no  abrupt  changes  in  the  pressure 
profiles  to  indicate  any  regions  of  discontinuity®  However, 
a  comparison  of  the  theoretical  steady  state  curve  for  a 
homogeneous  system  with  the  steady  state  curve  obtained 
indicates  a  substantial  discrepancy.  This  is  attributed  to 
the  lower  permeabilities  experienced  in  the  Berea  sandstone 
near  the  producing  boundary.  Inverting  the  Series  core  pro¬ 
vided  drastic  changes  in  the  pressure  profiles.  It  now  be¬ 
came  apparent  that  a  discontinuity  did  exist  somewhere  in 
the  core,  as  can  be  seen  in  Figures  26  and  27  of  Appendix 
IV,  However,  it  was  not  possible  to  determine  the  exact 
location  of  the  heterogeneity  from  the  graphical  results. 

A  further  inspection  of  the  results  for  the  Series  core  in¬ 
dicated  that  the  presence  of  the  Berea  fracture  did  not  notice 
ably  influence  the  shape  of  the  curves  obtained  in  this  region 


49 


As  a  result  it  was  concluded  that  the  Berea  sandstone- Alundum 
interface  had  no  effect  on  the  transient  results  obtained, 
and  the  discontinuity  detected  in  this  region  was  due  entirely 
to  the  difference  in  permeabilities.  With  the  exception 
perhaps  of  the  Indiana  limestone,  no  sharp  or  distinct  breaks 
in  the  pressure  profiles  were  discernible  for  the  results  ob¬ 
tained.  This  has  been  ascribed  to  the  fact  that  values  for 
the  pressure  are  known  only  for  specified  points  within  the 
core,  and  drawing  a  curve  through  these  points  would  tend  to 
average  or  smooth  out  the  effects  of  any  discontinuities  lo¬ 
cated  between  pressure  taps,  if  such  discontinuities  actually 
did  exist.  It  is  suggested  that  one  way  to  overcome  this  dif¬ 
ficulty  would  be  to  position  the  pressure  taps  closer  to¬ 
gether  in  a  region  of  known  heterogeneity.  One  pressure  tap 
would  be  located  in  the  high  permeability  region  while  the 
second  would  be  positioned  a  small  distance  away  in  the  region 
of  lower  permeability.  However,  care  should  be  taken  to  in¬ 
sure  that  the  presence  of  the  pressure  taps  would  not  markedly 
affect  the  transient  response.  Series  cores  with  sections 
having  greater  differences  in  permeabilities  than  those 
studied  here  would  provide  the  best  results. 

A  comparison  of  Figures  25  and  27  (Appendix  IV) 
indicates  that  for  the  sealed  boundary,  the  Series  core  ADS 
exhausts  itself  much  more  quickly  than  when  the  Alundum  sec¬ 
tion  is  located  in  the  upstream  position.  This  is  to  be 
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expected,  since  the  flow  from  the  high  permeability  Alundum 
for  the  Series  core  AUS  is  hindered  by  the  presence  of  a 
lower  permeability  Berea  sandstone  located  in  the  downstream 
position.  With  the  Alundum  section  in  the  downstream  posi¬ 
tion,  the  rate  of  flow  from  this  region  is  dependent  only  on 
the  back-pressure  applied  at  the  producing  face.  Although  the 
time  for  complete  exhaustion  should  theoretically  be  the 
same  for  both  cases,  any  specified  fraction  of  the  total 
available  gas  would  be  recovered  much  more  quickly  from  the 
Series  core  with  the  section  of  Alundum  in  the  downstream 
position.  Figures  25  and  30  (Appendix  IV)  serve  to  show  that 
the  sections  of  high  permeability  account  for  only  a  fraction 
of  the  total  pressure  drop  experienced  in  a  heterogeneous 
system.  Note  the  comparatively  small  pressure  gradients  in 
these  regions. 

Substantial  errors  are  suspected  in  the  values 
measured  for  the  variation  of  upstream  and  downstream  pres¬ 
sures  with  time.  The  crude  method  in  which  these  values 
were  obtained  leaves  much  to  be  desired,  especially  for  the 
higher  permeability  cores  where  changes  in  pressure  occurred 
much  more  rapidly.  This,  however,  is  a  limitation  of  the 
apparatus  used,  and  it  is  suggested  that  in  future  investi¬ 
gations  ,  a  measuring  device  be  installed  to  automatically 
sense  and  record  the  variations  of  both  the  upstream  and 
downstream  pressures  as  a  function  of  time.  The  gas  en¬ 
trapped  in  the  line  leading  from  the  pressure  tap  to  the 
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upstream  gauge,  and  the  free  gas  located  in  the  gap  created 
by  the  flange  and  core  surface  provided  a  second  source  of 
error,  since  this  gas  must  be  removed  through  the  core  in 
order  for  the  gauge  to  detect  a  change  in  pressure. 


Dimensionless  Quantities 


For  homogeneous  systems  the  theory  states  that  if 
the  parameters  influencing  transient  behavior  are  expressed 
in  dimensionless  form,  the  value  of  the  dimensionless  pres¬ 
sure  at  the  same  dimensionless  distance  and  dimensionless 
time  should  provide  a  single  and  unique  answer,  irrespective 
of  the  homogeneous  characteristics  of  the  core  and  the  nature 
of  the  gas  being  tested.  Figure  10  represents  an  attempt  to 
validate  this  concept.  A  dimensionless  time  of  3.0  was 
arbitrarily  chosen  and  the  values  of  real  time  for  each 
core  determined  from 
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where  <|>e  is  the  effective  gas  porosity  equal  to  <(>  (1-S  )  0 

The  values  of  viscosity  where  evaluated  at  P.  from  a  cor- 
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relation  available  in  the  literature 


(26) 


The  Visicorder 


charts  for  the  sealed  boundarv  case  were  then  read  at  these 
values  of  time  to  obtain  the  pressure  profiles,  and  hence 
the  dimensionless  pressure  distribution.  The  definition  of 
dimensionless  pressure  as  presented  here  is  valid  only  for  the 
constant  terminal  pressure  case.  Since  this  was  not  achieved 
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in  the  high  permeability  cores,  it  was  necessary  to  average 
the  downstream  pressure  for  the  time  interval  considered  in 
order  to  obtain  an  appropriate  value  for  p  .  Figure  11 
further  exemplifies  the  results  of  Figure  10.  In  this  in¬ 
stance,  however,  the  transient  behavior  for  various  times 
is  presented  at  the  same  position  in  each  core. 

The  cumulative  molar  production  n  from  each  core 

c 

at  any  time  t  was  evaluated  from  a  simple  material  balance 


nc  =  ni  -  n (t )  (40) 

where 

=  number  of  moles  initially  in  core 
n(t)  =  number  of  moles  remaining  in  core  at 
time  t 

Application  of  the  equation  of  state  for  a  gas  leads  to 
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where  subscript  "i"  refers  to  initial  conditions  anc 
script  "av"  refers  to  average  conditions  at  time  t. 
ing  the  results  of  equations  (23)  and  (41)  provides 
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where  P±  and  P0  are  equivalent.  Equation  (42)  was  then  used 
to  evaluate  the  dimensionless  cumulative  production  as  a 
function  of  dimensionless  time  for  each  core.  The  results 
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are  shown  in  Figure  12.  Note  that  although  P  and  P  re- 

av 

present  mean  pressures,  the  two  are  not  numerically  equal. 
Values  for  were  obtained  from  the  transient  pressure 
profiles  for  the  sealed  boundary  case.  By  definition  P 
was  then  evaluated  from 


2 

Compressibility  factors  were  obtained  from  data  available 
m  the  literature  for  nitrogen  gas(21).  Finally , correspond¬ 
ing  dimensionless  times  were  determined  from  equation  (39) 
by  replacing  P^^  with  P. 

From  the  theory  developed  for  homogeneous  systems, 
the  correlation  of  dimensionless  parameters  shown  in  Figures 
10  and  11  should  provide  a  single  curve.  An  inspection  of 
the  dimensionless  pressure  profile  for  the  Indiana  limestone 
in  Figure  10  shows  why  this  need  not  necessarily  be  true 
for  heterogeneous  porous  systems.  Although  the  dimension¬ 
less  correlation  has  incorporated  all  the  basic  variables 
influencing  the  transient  response,  no  parameter  has  been 
introduced  which  will  account  for  the  variation  of  perme- 
ability  from  point  to  point  within  each  core.  However,  it 
is  this  very  same  permeability  variation  which  is  respon¬ 
sible  for  the  unique  shape  of  the  pressure  profile  obtained 
for  the  limestone  core.  To  clarify  this  further,  consider 
the  results  obtained  for  the  Series  core  in  Figure  11.  The 
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point  of  interest  is  the  center  of  the  core,  which  is  exact¬ 
ly  the  same  regardless  of  the  position  of  the  Alundum  sec¬ 
tion.  A  comparison  of  the  results  for  the  Series  core  AUS 
and  ADS  indicates  a  large  discrepancy  in  the  values  of  dimen¬ 
sionless  pressure  for  the  same  dimensionless  time.  Since  the 
initial  pressures  were  very  nearly  the  same  for  both  tests, 
it  may  be  concluded  from  equation  (39)  that  the  real  times 
under  consideration  are  also  the  same.  However,  it  has  al¬ 
ready  been  shown  that  the  discrepancy  in  pressures  for  the 
Series  core  can  be  attributed  directly  to  the  difference  in 
orientation  of  the  high  permeability  section  with  respect  to 
the  producing  boundary.  Consequently,  before  the  transient 
response  of  a  heterogeneous  system  can  be  defined  uniquely, 
a  correlating  parameter  is  required  which  will  not  only  in¬ 
corporate  a  variation  in  permeability ,  but  will  also  account 
for  the  orientation  of  the  heterogeneities  with  respect  to 
a  fixed  reference. 

The  Series  Core  AUS  and  Berea  Sandstone  2  in  Figure 
10  provided  results  more  in  keeping  with  those  expected  for 
a  dimensionless  correlation.  The  dimensionless  pressure 
profiles  compare  favorably  within  experimental  error.  A 
comparison  of  Figures  25  and  28b  (Appendix  IV)  shows  that 
the  pressure  histories  are  very  similar  for  both  cores.  An 
inspection  of  Table  1  reveals  that  this  can  be  attributed  to 
the  similarity  of  the  permeability  distributions. 
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Figure  12  tests  the  applicability  of  Van  Everdingen 
and  Hurst's  dimensionless  cumulative  production  term  to  the 
experimental  results  obtained.  As  predicted  by  the  theory, 
values  appear  to  be  uniquely  defined  in  the  early  stages  of 
the  producing  life  of  the  cores.  Very  little  data  were  avail¬ 
able  to  substantiate  this  early  production  history  because 
of  the  small  values  of  real  time  involved,  especially  for 
the  high  permeability  cores.  The  presence  of  large  residual 
fluid  saturations  greatly  reduced  the  effective  porosity  of 
the  cores  to  gas.  As  a  result  the  cores  were  exhausted  more 
quickly  than  they  would  have  been  if  only  the  gas  phase  was 
present.  The  predicted  deviations  from  the  single  line  are 
due  to  the  influence  of  the  sealed  external  boundaries.  The 
irregular  shapes  of  the  curves  obtained  for  the  high  per¬ 
meability  cores  can  be  ascribed  to  the  fact  that  these  tests 
did  not  in  actuality  simulate  the  constant  terminal  pressure 
case,  a  prerequisite  for  the  use  of  this  correlation.  The 
results  for  the  tight  cores,  however,  show  more  of  a  resem¬ 
blance  to  the  curves  predicted  by  the  theory.  For  reasons 
previously  discussed,  there  is  a  large  discrepancy  in  the 
results  obtained  for  the  Series  core  AUS  and  ADS.  A  compari¬ 
son  of  the  two  curves,  however,  does  serve  to  indicate  the 
extent  of  the  error  that  may  be  introduced  by  applying  this 
theory  to  a  heterogeneous  system. 
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Radius  of  Drainage 

Results  from  equations  (28)  and  (29)  were  com¬ 
pared  with  values  for  the  radius  of  drainages  and  stabiliza¬ 
tion  times  determined  experimentally.  Knowing  the  values  of 
Pi  and  Pf'  ec2uation  (27)  was  used  to  calculate  the  value  for 
the  pressure  at  the  radius  of  drainage.  The  Visicorder  charts 
then  provided  a  means  for  establishing  the  time  required  for 
the  pressure  to  reach  this  level,  at  the  different  pressure 
taps  along  the  core.  Appropriate  values  of  Pf  for  the  high 
permeability  cores  were  again  obtained  by  averaging  the 
downstream  pressure  over  the  time  interval  considered.  The 
results  for  the  Indiana  limestone  are  shown  in  Figure  13. 
According  to  equation  (29)  a  log-log  plot  of  stabilization 
time  versus  radius  of  drainage  should  yield  a  straight  line 
with  a  slope  of  2,  if  y,  $ ,  k ,  and  P  are  constants.  These 
latter  restrictions  are  inherent  in  the  original  derivation 
of  equation  (29).  For  this  investigation,  however,  the 
above  limitations  were  mainly  overcome  by  using  the  effective 
permeability  of  only  the  section  of  core  through  which  the 
transient  had  traversed.  Equation  (38)  was  used  to  deter¬ 
mine  these  effective  permeabilities  as  the  transient  pro¬ 
gressed  from  point  to  point  within  each  core.  Curve  B  of 
Figure  13  was  determined  in  this  manner. 

A  further  inspection  of  equation  (29)  exemplifies 
the  fact  that  the  stabilization  time  is  directly  proportional 
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to  the  amount  of  gas  present,  represented  by  <j> ,  and  inversely 
proportional  to  the  rate  at  which  the  gas  may  be  removed, 
represented  by  k.  This  relationship  between  the  stabiliza¬ 
tion  time  and  the  porosity  and  permeability  is  tested  in 
Figure  14  assuming  a  constant  viscosity  and  mean  pressure. 
Variations  in  y  and  P  are  accounted  for  in  Figure  15. 
Theoretically  speaking,  all  values  obtained  for  Figure  15 
should  fall  somewhere  on  the  45°  line.  The  influence  of  per¬ 
meability  alone  on  the  stabilization  time  is  shown  in  Figure 
16.  As  expected,  the  stabilization  times  decrease  for  in¬ 
creasing  permeabilities,  all  other  factors  being  equal.  The 
permeabilities  obtained  for  Figure  16  were  determined  on  a 
section  by  section  basis  for  each  core  using  equation  (38) . 

The  stabilization  times  were  determined  from  equation  (29). 
Note  that  in  all  of  the  above  cases,  the  porosity  referred 
to  is  the  effective  porosity  to  gas. 

Figure  13  shows  that  a  large  discrepancy  exists 
between  the  theoretical  and  experimental  values  obtained 
for  the  stabilization  times.  This  is  believed  to  be  due  to 
the  relatively  large  reaction  time  of  the  experimental  ap¬ 
paratus  to  changes  in  pressure,  and  to  the  fact  that  it  is 
not  possible  to  obtain  an  instantaneous  pressure  drop  at  the 
producing  boundary.  Because  the  values  for  the  stabiliza¬ 
tion  time  are  so  small,  the  presence  of  the  initial  excess 
in  the  pressure  lines  leading  to  the  transducers  becomes 
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Figure  14  Influence  of  <j>/k  on  the  Radius  of  Drainage 
and  Stabilization  Time 
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Figure  15  Influence  of  y^/Pk  on  the  Radius  of  Drainage 
and  Stabilization  Time 
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Figure  16 


Influence  of  k  on  the  Radius  of  Drainage  and 
Stabilization  Time 
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significant.  This  is  especially  true  for  the  high  perme¬ 
ability  cores,  since  the  gas  entrapped  in  these  lines  must  be 
removed  through  the  core  before  a  pressure  differential  can 
be  detected  by  a  transducer.  Figure  17  is  an  example  of  a 
case  where  the  experimental  results  contradicted  those  pre¬ 
dicted  by  the  theory.  Note  that  the  relative  positions  of 
the  experimental  curves  are  the  direct  opposite  of  the  results 
expected  from  the  theory.  For  reasons  to  be  discussed  later, 
it  was  not  possible  to  obtain  a  value  for  the  stabilization 
time  at  all  the  pressure  tap  locations  in  each  core.  How¬ 
ever,  enough  data  were  available  to  show  that  although  the 
actual  values  for  the  stabilization  time  did  not  compare,  the 
general  shape  of  the  experimental  curve  was  accurately  pre¬ 
dicted  by  the  theory  (see  Figure  13) .  Consequently  it  was 
concluded  that  the  discrepancy  in  experimental  and  theoreti¬ 
cal  results  was  due  to  the  limitations  of  the  equipment  and 
not  the  theory. 

The  values  obtained  for  stabilization  times  at  the 
external  boundary  for  both  the  homogeneous  and  heterogeneous 
systems  are  the  same, since  this  calculation  is  based  on 
identical  effective  permeabilities.  Curve  B  in  Figure  13 
shows  that  when  the  effective  permeability  up  to  the  radius 
of  drainage  is  lower  than  the  overall  permeability  of  the 
core,  the  values  of  stabilization  time  for  the  heterogeneous 
system  fall  above  the  straight  line.  Values  below  the 
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straight  line  would  indicate  a  higher  effective  permeability 
up  to  the  radius  of  drainage,  and  a  point  of  intersection 
would  mean  that  the  permeabilities  are  the  same.  The  Series 
core  in  Figure  17  shows  that  in  one  instance  the  effective 
permeability  is  always  higher  than  the  overall  permeability, 
and  in  the  second  case  the  opposite  holds  true. 

The  derivation  of  equation  (29)  is  based  on  the 
fact  that  the  dimensionless  time  is  always  the  same  at  the 
radius  of  drainage  for  a  homogeneous  system.  This  concept, 
however,  was  not  violated  in  the  application  of  equation  (29) 
to  a  heterogeneous  system.  In  effect, a  number  of  homogeneous 
cores  were  assumed  having  the  same  permeability  as  the 
effective  permeabilities  determined  from  equation  (38) .  The 
stabilization  times  were  then  evaluated  for  these  hypotheti¬ 
cal  cores  at  the  desired  locations,  and  the  results  compiled 
to  obtain  a  history  for  the  heterogeneous  system  being 
analyzed.  A  plot  of  the  effective  permeabilities  used  and 
the  corresponding  stabilization  times  evaluated  is  shown  in 
Figure  16.  Note  that  all  the  points  fall  on  the  predicted 
straight  line  having  a  slope  of  minus  one. 

An  attempt  to  correlate  the  factors  influencing 
transient  behavior  with  the  experimentally  measured  stabili- 
zation  times  is  shown  in  Figures  14  and  15.  The  experimental 
data  obtained  in  either  case  do  not  even  remotely  resemble 
the  straight  line  predicted  by  the  theory.  This,  however, 
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is  due  to  the  unreliability  of  the  data  and  not  the  theory, 
especially  in  the  case  of  the  high  permeability  cores .  The 
underscored  data  points  in  Figures  14  and  15  represent  re¬ 
sults  for  the  same  radius  of  drainage  under  identical  cir¬ 
cumstances  .  Note  the  large  discrepancies  in  the  values  of 
stabilization  times  obtained,,  As  can  be  seen,  the  tight 
cores  provided  the  higher  stabilization  times  and  the  least 
amount  of  scattering  of  data  points.  Figure  14  is  valid  only 
for  the  ratio  of  viscosity  to  mean  pressure  used  in  this  in¬ 
vestigation,  The  influence  of  these  two  variables  on  the 
stabilization  time  is  included  in  Figure  15 „  However,  since 
the  values  of  viscosity  and  mean  pressure  were  nearly  the 
same  for  all  tests  and  since  the  variation  throughout  the 
course  of  a  run  was  small,  the  two  plots  are  very  similar. 

Numerical  Solution 

A  solution  to  equation  (32)  was  effected  in  terms 
of  pressure  squared  in  an  attempt  to  determine  the  influence 
of  assuming  an  average  pressure,  P,  on  the  final  results, 
p  in  this  instance  was  defined  by  the  finite  difference  ex¬ 
pression 

P  =  —  (pr+]_  s  +  i  +  Pr,s+1  +  Pr-l,s+l  +  Pr+l,s 

6 

+  P  ,  ) 

r , s  r- 1 , s 


+  P 


(44) 


♦ 


where  the  asterisk  denotes  values  from  the  previous  iteration. 
Because  of  limitations  in  the  experimental  apparatus  and  the 
numerical  solution,  the  influence  of  using  P  was  not  accurately 
determined b 

The  experimental  results  were  expressed  in  terms  of 
dimensionless  parameters  and  a  comparison  made  with  the  numeri¬ 
cal  solution.  Figures  18  through  21  inclusive  compare  the 
dimensionless  pressure  profiles  obtained  experimentally  with 
those  determined  numerically  for  the  same  dimensionless  times. 
The  solutions  for  the  heterogeneous  and  homogeneous  systems 
were  obtained  by  utilizing  the  finite  difference  approxima¬ 
tions  (11-10)  and  (11-11)  (Appendix  II) ,  respectively.  Equa¬ 
tion  (39)  was  used  to  evaluate  the  dimensionless  times. 

Figures  18  and  19  for  the  Indiana  limestone  simu¬ 
late  the  ideal  case  of  "Constant  Terminal  Pressure".  Be¬ 
cause  of  its  high  permeabilitv ,  the  Series  core,  on  the 
other  hand,  was  subjected  to  a  varying  pressure  at  the  produc¬ 
ing  boundary.  This  was  simulated  numerically  by  performing 
a  three-point  Lagrangian  interpolation  on  the  experimental 
data  obtained  for  the  variation  of  the  downstream  pressure 
with  time.  Consequently ,  as  shown  in  Figures  20  and  21,  the 
experimental  and  numerical  values  of  pressure  are  identical 
at  this  point. 

In  all  cases  the  length  of  the  transient  period 
predicted  by  the  numerical  solution  varied  significantly 
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Figure  13  Comparison  of  Experimental  and  Numerical  Results 


Indiana  Limestone  -  Sealed  External  Boundary 
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Figure  19  Comparison  of  Experimental  and  Numerical  Results 
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Dimensionless  Distance  x/L 

Figure  20  Comparison  of  Experimental  and  Numerical  Results 
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from  the  experimental  results  obtained.  As  can  be  seen  in 
Figures  18  and  20  for  the  case  of  a  constant  pressure  at  the 
external  boundary,  the  numerical  solution  converges  to  the 
steady  state  in  a  much  shorter  time  than  would  be  expected 
on  the  basis  of  the  experimental  results.  For  the  case  of 
the  sealed  external  boundary  shown  in  Figures  19  and  21,  the 
numerical  solution  predicts  a  shorter  time  required  to  effect 
complete  exhaustion  of  the  core.  As  an  example,  consider 
the  dimensionless  time  of  1.734  in  Figure  19.  Note  that  the 
average  dimensionless  pressure,  as  predicted  by  the  numerical 
solution,  is  approximately  0.35,  while  for  the  same  time  the 
experimental  data  provides  an  average  dimensionless  pressure 
of  approximately  0.60.  The  reasons  for  this  discrepancy  in 
the  values  of  time  is  believed  to  be  due  mainly  to  the  in¬ 
herent  limitations  of  the  experimental  apparatus .  For  the 
limestone  core  the  numerical  predictions  are  based  on  the 
boundary  condition  stating  that  an  instantaneous  pressure 
drop  is  applied  to  the  producing  end  at  time  zero.  It  is, 
of  course,  impossible  to  completely  simulate  this  condition 
in  practice.  The  numerical  solution  also  fails  to  account 
for  the  finite  quantities  of  gas  entrapped  in  the  pressure 
lines  leading  to  the  transducers. 

The  prime  source  of  experimental  error,  however, 
is  believed  to  be  due  to  the  machining  of  the  pressure  taps 
and  the  corresponding  damage  caused  to  the  permeability  in 
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this  region,  If  the  reduction  in  permeability  is  of  a  sig- 
nificant  nature,  it  is  feasible  that  substantial  pressure 
gradients  could  exist  across  this  localized  region  during 
the  transient  life  of  the  core.  Since  the  analysis  of  the 
data  assumes  that  the  pressure  in  the  transducer  lines  is 
characteristic  of  the  actual  pressure  in  the  core  at  any 
time,  large  errors  could  result  in  the  subsequent  interpre¬ 
tation  of  the  Visicorder  charts .  Some  evidence  does  exist 
to  substantiate  this  hypothesis,,  In  a  number  of  instances 
negative  galvonometer  deflections  were  detected  on  the  Visi¬ 
corder  charts „  This  would  only  occur  if  the  downstream  pres¬ 
sure  were  higher  than  the  upstream  pressure,  which  is  impos¬ 
sible,  However,  if  the  region  of  core  surrounding  the  down¬ 
stream  pressure  tap  was  damaged  and  the  corresponding  up¬ 
stream  tap  was  not,  it  is  possible  that  once  the  transient 
had  passed  through, a  larger  pressure  could  exist  in  the 
downstream  transducer  line  than  in  the  upstream  line.  This 
would  certainly  account  for  any  negative  deflections. 

Neglecting  the  solubility  of  nitrogen  in  kerosene 
could  also  lead  to  errors.  Thus,  it  is  recommended  that  in 
future  investigations  the  solubility  be  incorporated  in  the 
theory,  or  experimental  evidence  be  provided  to  substantiate 
neglecting  it  altogether. 

Although  no  agreement  could  be  obtained  on  the 
timing  of  the  profiles  shown  in  Figures  18  through  21,  a  re- 
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semblance  does  exist  in  the  general  shape  of  the  predicted 
and  experimental  curves.  The  correspondence  at  the  steady 
state  in  Figures  18  and  20  is  within  experimental  error  and 
for  the  case  of  the  Series  core  the  steady  state  profiles 
are  identical.  It  is  for  this  reason  that  the  discrepancy 
obtained  for  the  theoretical  and  experimental  transient  re¬ 
sults  was  attributed  to  the  apparatus  used  and  not  the  theory 
applied.  Note  that  no  attempts  were  made  to  account  for  the 
effects  of  slip  or  turbulence  on  the  theory  used  to  describe 
the  transient  behavior.  As  would  be  expected,  there  is  a 
time  correspondence  between  the  predicted  pressure  profiles 
obtained  for  the  homogeneous  system  with  those  obtained  for 
the  heterogeneous  core. 

Although  the  large  discrepancies  in  the  experi¬ 
mental  and  numerical  results  have  been  attributed  mainly  to 
the  equipment  used,  it  is  important  to  note  that  the  method 
of  expanding  equation  (32)  before  finite  differencing  could 
lead  to  stability  problems  in  the  numerical  solution.  As  a 
result  of  this  procedure,  the  matrix  of  coefficients  obtained 
is  no  longer  symmetrical.  Consequently  inherent  restrictions 
are  imposed  on  the  maximum  size  of  At  and  Ax  in  order  to 
insure  stability.  No  attempts  were  made  to  analyze  the  ef¬ 
fects  of  these  limitations  on  the  numerical  results  obtained. 
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Influence  of  the  Residual  Fluid  Saturation 

The  results  of  this  investigation  were  compared 
with  the  results  of  a  similar  investigation  conducted  by 
C.  Fortems  ,  in  an  attempt  to  determine  the  effects  of 
the  presence  of  the  residual  fluid  saturation.  A  compari¬ 
son  was  obtained  for  the  transient  pressure  history  at  the 
same  point  in  the  same  core  with  and  without  the  presence 
of  the  second  immobile  phase.  Both  experimental  and  numeri¬ 
cal  results  are  shown  for  the  Indiana  limestone  and  Series 
core  ADS  in  Figures  22  and  23  respectively. 

It  is  important  to  note  that  Fortems  did  obtain  a 
reasonable  agreement  between  his  numerical  and  experimental 
results  for  the  tighter  cores,  especially  in  the  later 
stages  of  the  transient  period  where  changes  in  pressure 
with  time  occur  less  rapidly  (refer  to  Figure  22) .  As  can 
be  seen  in  Figure  23,  however,  his  work  also  provided  large 
discrepancies  for  the  high  permeability  cores,  for  reasons 
previously  discussed.  The  fact  that  one  investigation  pro¬ 
vided  some  comparable  results  while  the  second  provided  none 
at  all  is  attributed  to  the  presence  of  the  residual  fluid 
saturation.  Since  the  fluid  cannot  be  completely  displaced 
from  the  finer  regions  of  the  porous  network,  a  reduction  in 
^ff^ctive  permeability  to  gas  is  experienced.  Conseguently 
it  is  feasible  that  the  kerosene  would  imbibe  into  the  tight 
regions  surrounding  the  machined  pressure  taps,  and  would 
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Figure  22  Influence  of  SQr  on  the  Transient  Behavior 
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Figure  23  Influence  of  S  on  Transient  Behavior 
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not  be  dispelled  in  the  subsequent  displacement  by  nitrogen 
gas.  As  a  result  the  reduction  in  permeability  in  these 
localized  regions  would  be  further  accentuated  by  the  pre¬ 
sence  of  the  liquid  phase. 

Even  more  startling  is  the  fact  that  no  agreement 
could  be  obtained  for  the  numerical  solutions  from  both  in¬ 
vestigations  on  a  dimensionless  basis.  Again  refer  to  Fi¬ 
gures  22  and  23.  Since  the  dimensionless  parameters  fail  to 
account  for  any  heterogeneities,  it  is  suggested  that  this 
may  be  a  reason  for  the  deviation.  As  can  be  seen  in  Fi¬ 
gure  19  from  the  numerical  solutions,  the  pressure  at  any 
point  is  more  dependent  on  the  surrounding  permeability  at 
this  point  than  on  the  overall  permeability  of  the  core.  For 
example, a  comparison  of  the  homogeneous  and  heterogeneous  so¬ 
lutions  at  a  dimensionless  distance  of  0.1667  and  dimension¬ 
less  time  of  0.867  yields  a  discrepancy  of  approximately 
25  percent  in  the  values  of  pressure  obtained.  However,  in 
both  cases  the  overall  effective  permeability  used  for  the 
core  was  nearly  the  same.  Consequently,  if  the  presence  of 
the  residual  fluid  saturation  were  to  reduce  the  permeability 
of  the  core  more  in  one  region  than  in  another  due  to  satur¬ 
ation  gradients  etc.,  the  overall  permeability  distribution 
would  be  altered,  and  a  discrepancy  in  numerical  results  at 
any  one  point  could  be  expected. 
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As  has  already  been  shown  in  Appendix  I ,  the  study 
of  gas  flow  in  the  presence  of  a  residual  fluid  saturation 
is  a  special  case  of  multiphase  flow.  Consequently,  a  more 
accurate  analysis  of  the  experimental  data  obtained  would 
include  the  concepts  of  total  compressibility  and  total  mobil¬ 
ity.  Since  only  the  gaseous  phase  was  flowing,  the  theory 
developed  for  single  phase  gas  flow  was  sufficient  to  describe 
the  actual  system  as  far  as  the  mobility  was  concerned.  How¬ 
ever,  a  number  of  approximations  were  required  in  order  to 
establish  the  total  compressibility  of  the  system  as  the  in¬ 
verse  of  the  pressure.  Among  these  were  the  assumptions  that 
the  solubility  of  nitrogen  in  kerosene  was  negligible,  and 
that  the  variation  of  the  gas  compressibility  factor  with 
pressure  was  small.  For  the  range  of  pressures  considered 
(0-65  psia) ,  the  approximations  were  valid  and  little  error 
was  introduced  by  the  application  of  single  phase  gas  flow 
theory.  Note,  however,  that  for  pressures  normally  en¬ 
countered  in  gas  reservoirs,  the  effects  of  gas  solubility, 
liquid  compressibility,  and  variations  in  the  gas  compres¬ 
sibility  factor  are  significant. 
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CONCLUSIONS  AND  RECOMMENDATIONS 

Conclusions 

1*  For  the  cores  tested  the  presence  of  the  residual  liquid 
saturation  increased  the  discrepancy  between  the  experi¬ 
mental  and  numerical  results  obtained. 

2.  The  unsteady  state  gas  flow  theory  appears  to  be  valid 
for  the  flow  of  gases  through  porous  media  in  the  pre¬ 
sence  of  a  residual  liquid  saturation  for  the  cores 
analyzed . 

3.  A  better  delineation  of  the  properties  and  some  modifica¬ 
tions  to  the  theory  are  required  to  account  for  the  dis¬ 
continuities  and  the  sequence  in  which  these  occur,  be¬ 
fore  a  heterogeneous  system  can  be  accurately  analyzed 

on  a  dimensionless  basis. 

4.  Serious  errors  can  result  when  applying  the  dimension¬ 
less  cumulative  production  term  to  a  heterogeneous  sys¬ 
tem  by  using  an  effective  overall  permeability.  This 
is  especially  true  when  the  variations  in  permeability 
from  point  to  point  within  the  porous  medium  are  large. 

5.  Although  the  stabilization  time  equation  was  derived 
for  a  homogeneous  system,  it  appears  to  be  applicable 
to  heterogeneous  porous  media  in  the  presence  of  a 
residual  liquid  saturation,  if  the  effective  sectional 
permeability  up  to  the  radius  of  drainage  is  used  in  the 


82 


calculations.  The  practical  applications  of  this  equation, 
however,  are  limited  by  the  restrictions  imposed  on  it 
through  the  boundary  conditions.  Note  that  equation  (28) 
can  no  longer  be  used  for  heterogeneous  systems,  since 
the  radius  of  drainage  must  be  known  in  order  to  obtain 
an  appropriate  value  for  the  permeability. 

6.  The  numerical  techniques  developed  for  describing  the  tran¬ 
sient  flow  of  gas  through  a  heterogeneous  porous  system 
appear  to  be  valid  for  practical  situations.  Although 
some  error  may  be  due  to  the  stability  problems  encoun¬ 
tered  m  the  numerical  solution,  the  large  discrepancies 
obtained  in  the  experimental  and  numerical  results  may  be 
attributed  mainly  to  the  limitations  of  the  equipment  used. 

7.  The  internal  fracture  or  core  interface  introduced  by 
simulating  two  different  permeabilities  in  series  did  not 
noticeably  affect  the  transient  response 0 

8.  The  mechanical  end  effects  were  eliminated  by  shearing 
the  core  faces  as  opposed  to  cutting.  All  discontinuities 
detected  were  thus  attributed  to  natural  heterogeneities 
or  damage  by  other  means. 

Recommendations 

1.  Only  long  and  tight  cores  should  be  used  in  an  experimental 
transient  analysis,  since  these  provide  the  best  results. 

2.  Automatic  measuring  devices  should  be  installed  at  both 
ends  of  the  core  in  order  to  obtain  an  accurate  and  perma¬ 
nent  record  of  the  pressure  variation  with  time  at  these 
locations . 
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3.  Attempts  should  be  made  to  eliminate  the  effects  of 
machining  the  pressure  taps. 

4.  The  effects  of  heterogeneities  may  be  more  accurately 
determined  if  the  pressure  taps  are  located  closer  to¬ 
gether  in  regions  of  known  discontinuity.  Also  it  is 
suggested  that  future  series  cores  be  simulated  with 
sections  having  a  greater  variation  in  permeability 
than  the  one  analyzed  here. 

5.  The  lengths  of  lines  leading  to  pressure  sensing  devices 
should  be  reduced  to  a  minimum.  Attempts  might  also  be 
made  to  account  for  the  gas  entrapped  in  the  void  space 
created  by  the  flange  and  core  interface. 

6.  Future  analyses  might  also  incorporate  the  effects  of 
the  solubility  of  gas  in  the  immobile  liquid  phase. 

7.  The  stability  problems  encountered  in  this  investigation 
concerning  the  numerical  solution  for  equation  (32)  can 
be  overcome  if  the  matrix  of  coefficients  obtained  from 
the  finite  difference  approximations  is  symmetrical. 

This  can  be  achieved  if  the  diffusivity  equation  for  gas 
flow  is  finite  differenced  in  the  form  expressed  by 
equation  (32)  .  It  is  recommended  that  this  procedure 

be  followed  in  future  investigations. 

8.  No  attempts  were  made  to  account  for  the  porosity  varia¬ 
tions  and  liquid  saturation  gradients  within  the  cores. 
Consequently  future  investigations  might  consider  varia¬ 
tions  in  <p  and  S 

or 
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NOMENCLATURE 

Where  possible,  standard  letter  symbols  have  been 
used  as  outlined  by  the  Society  of  Petroleum  Engineers  ( 41 } . 
Any  deviations  from  this  set  standard  have  been  duly  noted 
throughout  the  text  of  this  thesis. 
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A  P  P  E  N  D  I  X  I 


THE  CONCEPT  OF  MULTIPHASE  FLOW 
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The  equations  governing  the  simultaneous  flow  of 
gas,  oil,  and  water,  neglecting  the  effects  of  gravity  and 
capillary  pressure,  are  respectively 
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In  those  cases  where  the  pressure  and  saturation  gradients 

are  small,  the  vector  products  VP-VP,  VP -VS  ,  and  VP • VS  are 

o  w 

small  compared  to  VP,  vSo,  and  VSw.  For  such  instances 
( 31 ) 

Martin  has  shown  that  equations  (1-1)  through  (1-3)  can 

be  combined  into  the  single  expression 
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In  effect  c^_  represents  the  total  compressibility  of  the  fluid 
and  A  the  total  mobility. 

However ,  for  the  problem  being  considered  there  is 
no  water  phase  present  and  only  the  gas  phase  is  mobile. 
Equation  (1-6)  consequently  reduces  to 


k 

=  _a 


(1-7) 


and  if  the  solubility  of  the  gas  in  oil  is  considered  negli¬ 
gible*  , equation  (1-5)  reduces  to 


9B  S  9B 
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Applying  the  basic  definition  of  compressibility 

1  9V 

c  =  -  “  (— )T 
V  9P  1 

equation  (1-8) can  be  further  simplified  to 


(1-8) 


(1-9) 


c,  =  S  c  +  S  c 
t  O  o  g  g 
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*  To  date  no  experimental  data  exist  relating  the  solubility 
of  nitrogen  gas  in  kerosene  as  a  function  of  temperature 
and  pressure.  Consequently  it  is  difficult  to  draw  any 
conclusions  as  to  the  extent  of  the  error  introduced  by 
making  this . assumption .  An  attempt  was  made  to  calculate 
the  solubility  of  air  in  kerosene  to  obtain  some  form  of  a 
criterion  for  nitrogen.  The  results  showed  that  the  volume 
of  gas  dissolved  in  the  kerosene  could  be  as  high  as  18%  of 
the  free  gas  occupying  the  remaining  pore  space.  However, 
the  method  of  analysis and  the  data  employed  are  of  a 
dubious  nature,  and  consequently  no  concrete  conclusions 
could  be  drawn . 
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It  can  also  be  shown  with  the  aid  of  equation  (1-9)  that  for 
the  range  of  pressures  considered  in  this  investigation ,  c 
is  equal  to  the  inverse  of  the  pressure.  Consider  the  follow¬ 
ing  set  of  typical  values: 


Therefore 


S 

g 


0.55 


c 

o 


1.42  x  10  2  psi 
1  1 


0.45  c  =  -  =  _ 

^  P  65  psia 

=  1.54  x  10~2  psi"1 

=  0.55(1.42  x  10“5)  +  0.45(1.54  x  10"2) 

0.00781  x  10  2  +  6.93  x  10  2  psi  1 


* 


Since  the  compressibility  of  the  oil  is  negligible  compared 
to  the  compressibility  of  the  gas,  the  total  compressibility 
of  the  fluid  can  be  approximated  by  equation  (1-11)  intro¬ 
ducing  very  little  error 


c 


t 


(1-11) 


At  this  point  it  should  also  be  mentioned  that  the  compres¬ 
sibility  of  the  porous  system  has  been  neglected  in  equation 


This  value  for  the  compressibility  of  kerosene  was  obtained 
by  the  method  of  Trube [  >  ,  assuming  a  pseudo-reduced  pres¬ 

sure  of  1.00  in  order  to  obtain  an  approximate  value  for  the 
pseudo-reduced  compressibility.  Consequently  the  final  re- 
sul"t  itself  is  an  approximation.  It  does,  however,  serve 
to  illustrate  the  range  of  values  to  be  expected. 
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(15).  A  typical  value  of  c^  for  a  porosity  of  20  percent^^) 

*  ""6  "-I 

is  3.7  x  10  psi  .  A  comparison  of  this  value  with  those 
previously  mentioned  indicates  that  the  compressibility  of 

the  porous  material  adds  very  little  to  the  net  compressibility 
of  the  system. 

Substituting  equations  (1-7)  and  (1-11)  into  equa¬ 
tion  (1-4)  and  replacing  S  by  (1  -  S  )  provides 

g  o 
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- _ g _ 


k  P 
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9t 


(1-12) 


Now  multiply  both  sides  of  equation  (1-12)  by  2P, 
assuming  an  average  pressure  P  for  the  inverse  of  the  com¬ 
pressibility,  and  replacing  <j>  (1-SQ)  by  the  effective  poro¬ 
sity  to  the  gas  phase,  equation  (1-12)  becomes  (the  subscript 
g  is  understood) 


2  4*  y  9  P 

2PV  P  =  —  2P  —  (I-13) 

kp  at 

Consider  the  mathematical  relationship 

2  2  9 

VP  =  2 VP • VP  +  2PV  P  (I-14a) 

If  now  it  is  assumed  in  equation  (I-14a)  that  VP -VP  is  neg¬ 
ligible,  as  was  done  in  the  derivation  of  equation  (1-4) ,  the 
result  is 

2  2  2 
V  Pz  *  2P V  p 


(I-14b) 
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Substituting  this  relationship  into  equation  (1-13)  and 
simplifying  provides 

22  8p2 

VP  =  —  -  (1-15) 

kP  at 

A  comparison  of  this  equation  with  equation  (22)  of  the 
Theory  for  single  phase  gas  flow  shows  that  the  two  are 
identical . 
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appendix 

I  I 

DERIVATION  OF  THE  FINITE 

DIFFERENCE 

APPROXIMATION 


II-  2 


Defining 


W(P)  = 


M  (P) 


(II-l) 


and  applying  the  mathematical  equivalent. 
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expansion  of  equation  (32)  of  the  Theory  provides 


W(P)  k  (x) 


2  2 
9  P 


9x 


—  +  k  (x) 


3P2  9P  2  dW(P) 


9x  9x  dP‘ 


(II-2) 
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(II-3) 


The  underscored  terms  are  treated  as  averages  in  the  implicit 
computations . 

Substituting  the  following  finite  difference  approxi¬ 
mations  into  equation  (II-3) 
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the  following  definitions 
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provides 
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It  should  be  noted  that  equations  (II-4b) ,  (II-4c)  and  (II-4e) 

represent  the  implicit  Crank-Nicholson ( 29 ^  finite  difference 
approximations  centered  at  time  s+^.  The  running  variable 
for  position  is  represented  by  r  and  the  asterisk  denotes 
the  values  of  P  at  time  s+1  as  determined  from  the  previous 
iteration.  In  equation  (II-6) 
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If  equation  (II-6) 

is 

now  multiplied  by  2 (Ax) 

2/k (x) W (P )  (A) 

and  equations  (II-8a)  and  (II-8b)  utilized,  the  result  is 
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Expanding  the  subscripted  A's  and  collecting  like  terms  pro- 
vides  the  final  result 
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If  the  viscosity  and  permeability  are  constant,  equation 
(1I_10)  reduces  to  the  familiar  finite  difference  approxima¬ 
tion  for  a  linear  parabolic  partial  differential  equation. 
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The  method  of  solution  now  lends  itself  to  evaluat¬ 
ing  equation  (11-10)  for  every  point  within  the  grid  spacing, 

and  solving  the  resulting  set  of  linear  equations  simultan- 

2 

eously.  The  values  for  P  at  time  s+1  represent  the  un- 

2 

All  the  values  for  P  at  time  s  are  known  and  all 


knowns . 
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the  coefficients  in  equation  (11-10)  can  be  evaluated.  The 
values  of  D(P)  were  computed  directly,  as  an  analytical  ex¬ 
pression  for  viscosity  as  a  function  of  pressure  was  avail¬ 
able  from  the  literature ^26 ^ . 

Special  forms  of  equation  (11-10)  are  required  at 
both  ends  of  the  grid  system,  depending  upon  the  boundary 
conditions  applied.  For  cases  where  the  pressure  at  the 
boundary  is  explicitly  defined,  either  P2  (downstream 

^0  1  ^  St  1 

2 

boundary)  or  pr+1,s+1  (upstream  boundary)  represent  known 
values.  At  a  sealed  external  boundary 


8x 


To  simulate  equation  (11-12)  numerically,  a  fictitious  boundary 
is  assumed  a  distance  Ax  past  the  actual  boundary.  If  r  re¬ 
presents  the  position  of  the  actual  sealed  boundary,  then  the 
conditions  which  prevail  at  r+l  are  identical  to  those  at  r- 1. 
More  specifically, 

r+l ,  s+1  =  Pr-l,s+l  (II-13a) 

and 

r+l ,  s  =  Pr-L,s  (II-13b) 

The  final  solution  to  the  resulting  tri-diagonal 
matrix  was  achieved  by  the  method  of  Thomas ^29^. 
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APPEND  IX  III 


COMPUTER  PROGRAM 
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The  final  solution  to  the  problem  outlined  in 
Appendix  II  was  achieved  by  means  of  a  computer  program, 
written  in  Fortran  IV  computer  programming  language,  and 
executed  on  the  IBM  OS/360/67  model  computer.  The  program 
is  designed  to  solve  equation  (32)  of  the  Theory  for  essen¬ 
tially  four  separate  boundary  conditions: 

1)  constant  terminal  pressures  at  both  boundaries; 

2)  the  constant  terminal  pressure  case  with  a  sealed 
external  boundary; 

3)  a  varying  sandface  pressure  with  a  constant  terminal 
pressure  at  the  external  boundary;  and 

4)  a  varying  sandface  pressure  with  a  sealed  external 
boundary . 

The  initial  condition  in  all  cases  was 

P(x,0)  =  Pi  (III-l) 

Although  written  primarily  to  accommodate  variations  in  vis¬ 
cosity  and  permeability,  the  program  is  flexible  enough  to 
solve  equation  (32)  for  a  homogeneous  system  and  a  constant 
viscosity,  requiring  only  minor  adjustments. 

The  data  required  is  as  follows: 

1)  initial  pressure  (psia) 

2)  residual  fluid  saturation  (fraction) 

3)  porosity  (fraction 


I  I 
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4)  flowing  sandface  pressure  (psia) * 

5)  atmospheric  pressure  (psia) 

6)  temperature  (°F) 

7)  total  length  (inches) 

8)  incremental  length  (inches) 

9)  total  time  (minutes) 

10)  incremental  time  (minutes) 

11)  permeability  distribution  as  a  function  of  dimensionless 
distance,  x/L  (md)** 

12)  the  variation  of  the  sandface  pressure  (psig)  as  a 
function  of  time  (minutes)** 

13)  the  criterion  for  convergence  (psi) 

The  choice  of  boundary  conditions  was  effected  by  means  of  a 

positive  (pressure  varies)  -  negative  (pressure  is  constant) 
coding  system. 

The  final  computer  output  consists  of  a  pressure 
and  dimensionless  pressure  (P/P^  distribution  for  each  in¬ 
cremental  time  step  At.  Also  included  are  several  statements 
describing  the  nature  of  the  problem  and  the  boundary  condi¬ 
tions  being  applied. 


cases  the  constant  flowing  sandface  pressure, 
where  required,  and  atmospheric  pressure  were  the  same. 

These  data  are  optional, and  depend  upon  the  problem 
being  solved  and  the  boundary  conditions  applied.  In 
cases  where  either  y  or  k  are  considered  constant,  it 
is  necessary  to  read  in  the  respective  values  separately. 


- 


Ill-  4 


Wherever  possible  standard ^9)  mnemonic  computer 

symbols  have  been  used.  Most  deviations  from  this  standard 

are  self-explanatory  within  the  program  itself.  Exceptions 

to  the  above,  however,  are  listed  below  for  further  clarifi- 
cation . 


PRSXP2 

PRSREQ 

PRSCS 

DELK2 

DELK1 

PRMGVT 

B  (M,N) 

D  (M) 

SUBRTC 

BCC 

PECHEK 

VCHECK 

TCHECK 

FCHECK 


Pf  at  time  t 
Pf  constant 


k (x/L) 

coefficients  of  the  tri— diagonal  matrix 
corresponding  vector  for  B(M,N) 

storage  locations  for  the  positive-negative 
codes 

codes  which  determine  the  execution  of 
certain  operations  depending 
on  their  respective  values 


It  should  be  noted  that  the  computer  program  as  presented 
here  represents  a  working  solution  to  the  problems  dis¬ 
cussed.  Absolutely  no  attempt  has  been  made  to  optimize 
the  calculations  by  available 


programming  techniques. 


CGMMCN/MANE/PRSXP2(3,105)  »  B  (  105,3)  » D  (  100)  ,  K 
CGMMGN/CALC8/5UBRTC(3  ) ,PRMG , BETA , bETAI ,R 

- I - CflMMQN/VlSC  S  ¥  /R  RSVI  £ «  TFM  ,  VCnECK, VISG 

COMMON /F  LISBON/ DEL  P  ,WOFPRS »OELPKG,F CHECK 
COMMON/ PERM/ PECHEK  »  DELLTH , DEL PRM , N 
- QflMFCN/PCLATE/PRSREu. TIM. TCHECK _ 

COMMON  PRS  t  3  » 1 0 5 ) 

INTEGER  R 

_ JB.I.ME.NS1U N  BCC  (3  )  ,PRSQ  (  100  ) 

REAL  PRS I tSATLRf PORT , PR SCS , D ELT I M , D E LLTHfTI MT , PORE, 7  EM 
1  »  L  Th , DEL  PRS 

- REALM  5. 1C)  (SUfikTCt  I  ) « 1= 1 «  3  ) _ _ 

1C  FORMAT ( IX, 3( 1X,F4.1 ) ) 

WRITE (6, 10) (SUBRTC(I) ,1=1,3) 

R  E  A  U  (  3  »  1 1  )  PRSI .SAT  LR .PORT . PR SOS  .TPM ,LTH 

11  FORM AT I IX * F7. 1 , 2 { IX , F  7 . 4 ) , 1 X , F7 . 1 , 1 X , F5« 1 , F7 . 3) 

WRITE  16, 11)  PRSIfSATLR,PCRT,PKSCSt7EM,LTH 
- READ ( 3 « 12 )  PELT  IM.Urt 1  TH. Ti MT _ 

12  FORMAT  (  1X,3(  iX,F6.3) ) 

WRITE (6, 12)  DELT  IM, DELL  TH , T I MT 
RcAu (5,13)  iaCC (1) ,1=1,2) 

13  FORMAT ( TX*2C 1X,F4.1 ) ) 

WRITE (6, 13)  (BCC<I),I-1,2) 

- RFAlM  5f  ?h  )  DFi  PRS _ _ _ 

2b  FORMAT (IX, F7. 4) 

WRITE(6,26)  DELPRS 

_ M= I lMl/DcL  LIM 

N  =  L  Ih/CELLTH 
T  I M  =  0  •  0 

_ EOKE=PCRT*( 1 .O-SATl  R) _ 

N  P  2  =  (\  +  2 
DO  30  k=  1 ,  NP2 

_ _PRS1 1,  R)  =  PRSI  _ 

PRS (2, R)  =  1 .0 
PRS ( 3  ,R)  =  PRS  I 

_ RE-SXP2  (  l.k)=PRSI**2.Q _ 

PRSXP2 ( 2  ,  R  )  =  1  •  0 

30  PRSXP2(3,k)=PRSI**2.G 

_ _ JBLfcXAJ  =  Q E L T I M  /  (  PCKE*DELLTH**2 .0 )  *0,000  632  8  3  _ 

IF(bCC(2) .GT.O.O)  GO  TO  111 
WR1TE(6, 23)  PRSI 

-23  FORMAT (IX, «TFE  PRESSURE  aT  THE  EXTERNAL  BOUNDARY  IS  On*. 
1 ' NS T ANT  AT  »  , F 7 . i , 1 P S I A « ) 

GO  TO  112 

111  WRITE(6,24) . . . 

24  FORMAT ( IX , ' THE  POROUS  MEDIA  IS  FINITE  WITH  A  SEALED  EX«, 
1  *  T ERN AL  BOUNDARY' ) 

112  IF(SUBRTC( 1) «LT.Q.C)  GO  TO  100 _ 

w  R  I T  E  (  6 , 1 4  ) 

14  FORMAT ( IX, • Tht  DOWNSTREAM  PRESSURE  VARIES  WITH  TIME') 

_ GO  TO  102  _ _ _ _ _ _ _ _ _ _  __ 

100  WRI’TE ( 6»  15)  PRi.es 

15  FORMAT ( IX , 1 THfc  DOWNSTREAM  PRESSURE  IS  CONSTANT  AT  *,F 

1  7  .  1 ,  «  P  S  I  A  '  )  
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102  If-  (  $CBRTC(2)  .01,0. C)  GC  TG  103 
RE AD ( b  »  16  )  PRMG 

_ 1.4  E.QRMAI  1 1 X » F  8 . 2 ) . . . 

^ll[(fcf17J  PKMG 

17  FORMAT ( IX THE  PERMEABILITY  CP  1  HE  POROUS  MEDIA  IS  CON*, 

_ 1..*  S  IAN  I  AT  *  .  F  8  .  /  .  *  M  i  i  1  ) _ _ 

GO  TO  104 
1C  3  WRITE(6,18) 

18  FORMAT  UX  ,  '  THE . PERMEABILITY  Ui  IiiE_jP.yRQ.US  MEDIA  VAR1  ES  *  , 

1*  hi  T  H  POSIT  ION*  ) 

1C4  IF(SUBRTCi3) .GT.Q.O)  GO  TO  10b 

_ READ (b .19)  ViSG _ _ 

19  P0RMA1  (  IX  ,  P  8  •  5 ) 

VvRITt{6,20)  VISG 

2JQL FORMAT (IX, *  I.UE . GAS  VISCOSITY . IS  CONSTANT  AT  UF6.5, *C* , 

1  *PS  *  ) 

GO  TO  106 

■1C  5  WRlTH6t2  1) _ _ _ _ _ ; 

21  FORMAT UX, ‘THE  GAS  VISCOSITY  IS  TEMPERATURE  AND  PRESSU', 
l'Rt  DEPENEENT' ) 

— — ™JP£C  HEK  =  0  •  0 . . 

VChlCk=O.C 
TCH ECK=0  •  0 

106  CO  200  K.  S  =  1  f  M _ _ _ _ 

T  I  M=  T  I  M+DtLI  IM 

AGO  IP ( 6CCI 1) .LT .0.0)  GG  TL  107 

_ CALL,  TERPCL _ _ _ _ _ _ _ 

F*SXP2(2,l)-= ( PRSkEO+PRSCS ) **2.0 
GC  TO  108 

107  P  R  S  X  P  2  4  2  .  i  )  =  PkSr.S**7-C _ _ _ _ 

1 C  8  R  =  2 

PRSXP2I  1»  1)  =  PRSXP2(2,  1) 

_ CALJ _ SUdELA _ _ _ _ _ _ _ : . . 

^yNSTP=l  .  +  CLlP 
C0NSTM=1 .-DELP 

_ bMt2)=—  2.CU1.  +  1./8E1A) _ _ _ 

P( 1,3J=CCNSTP 

Di 1 )=-PRSXP2( 1, 3) *CGNSTP+2  .*( 1.- 1 . /BETA ) *PRSXP2 ( 1 ,2)-P 
_ lR^£2i.XiJJli«lCG S  T M-P.RSXP2.L2i  1 1  *C0NSTM-0E1.PRM 

1P(  BCCU  )  .LT.O.O)  GO  TG  109 
R  =  N  +  1 

_ PftS(liR»l)  =  PKSU.R-l) _ _ 

PRS(3»R+1)=PRS( 3,R-1) 

PRSXP2I i,k+l)=PRSXP2( 1 »  R—  1  > 

_ P&SAP2 ( 3 . 8  tl )  =  PR  S  X  P  2 ( 3 . R- 1 )  _ _ 

CALL  SU8ETA 
BIN, 11=2. 

_ d  i  N »  2 )  =~2  .-*U.  +  1./6.LLAJ _ 

D(N) =-2.*PRSXP2 ( 1»  R-l )+2.*( 1 .-1. /BETA  )*PRSXP2 U  ,R  ) 

K  =  N 

GO  TC  110 _ 

1C9  R=N 

PRSXP2I2  ,NU)=PRSl**2 .0 
CALL  SG3ETA 


i  j  .  t  .  i  «.  .  I  V.  )  j  I  >'  J jc  )  1 1  S  j  1 
o  -M  I  ,u  »  ')  jAj-l 
(  ' ,  .  )  I  ,  v  .  /  i  V  .  ,i  _  i  j  1 
t  i  N  i  »  Mil  v 
'  • 

it'  1  .  h  I  £  '  i 

( i ,  ;  -j  1 1  h  n  t  oi 


i  I  I  I  1  ,  y  .  )  *  ^  I  1 

I  1  '  v  1  i  i  I  I  i  *  I 

j  ;  .  .  I  U  .  (  c  ;  I  -  J JC  )  ii  *J1 

UC  1  Id.  )  UM  ‘  *i 


M  .  f 


(  .  * ,  i  )  i  1  H  u 

tu:  i  .  I  .  »  )  J  I  x  *  ■* 

11  ,  v.  >CV  .  ’  I  i  f  t  A  i  )  f  I  I  0^ 

l  '  <H  '  1 
did  J  f  lJD 
t  j  )  j  I  1  •>  j  Mi 


ci  y.IcOJ<ilV  iT'tXii  f  AMHOHI  15* 

(Mi  J  j  .'M 

:  >  • ( - aj H J  j 4 
0  .  _  */M  . :HJV 
0.0*  A  JJHJT 

t  I  -  «  ■  ^  v  1  L 


i  1  J,  u-M  l  l  *M1  l 

ii  j  I  l  .  .  i  *  I  1  )  )  d  U  0  »> 

j  1H1  j  i  j J4J  d 

•  )  t  i  ,  /  ,  <  ^  t 

b  v>  I  J  I  JO 

••  .>o  i  -.1  i  ,  .  /  vc  i  LI 


601 

Mt  J  'At..  -  (  1  t  l  I  .  M  A  <i  *  i 

J  J  J*  .  i  -  i  l  d.Lj 

*  J  ,  -  .  1  .  C  JJ 

i 

<  1 d  JJ-  ( t « 1 )  i 


,  .•  l  I  \  .  - .  i  >  •  .  j  ♦  4  I  I  l  il  h  iAC  /  H-s{lK 

i  i  ,  r  -  it  /  u  » i  / 

J  (I.v-.TJ.  liLddHI 

l+A-H 

i 1 ~  •  t i  / t  •  •  I  tl)C  11 

i  :  *  »  >  v  -M  »  fc  ) 

l  1  -  t  ■  1  A  <  4  -  (  i  ♦  i,  i  >  .  1/  <J  M 
M  -  ,  J  M  A<  I  .  A.,  d 

i I t  H)  i 


I  I  \ .  1  * .  i  1  ' .  -  -  (  »  1  J 

01J  JT  JJ 


t>01 

j  .  id  i*i  (!♦/  tiJ'i'i'd  '( 


R  =  N- 1 

CGf\iSTP=l  .+DELP 

_  CQNSIM=1 «- C  L  l  P 

6  (H,  1  )=CGNSI  M 
B(R,2)=—  2.*(  1  .  +  1  ./BE™  ) 

- DIB.)  =-PRSXP^  (  1,R+2)*CCNSTP+2.*(1.-1./6ETA)*PRSXP2(  l.Kt 

1 1  )— PRSXP2 I  1 »  k  ) *CGKS  TM-PRSXP2 {  2,R  +  2) *CQNSTP-DELPRM 
K  =  N-  1 

-HQ  00  199  K=3,K  _ [ _ _ _ 

CALL  SUBETA 
CCIN;  S  TP  =  i  •  +  GE  lP 

_ CUKSTM=l.-DbLP _ _ _ 

B(R-1,1)=CGRSTM 
B(R-l»2)=-2.*( l.+l./BETA) 

aiR=1^3LA.«CONSTP  _ _ __ _ _ _ _ _ _ 

IS 9  L  (  R- 1 I PR SXP2  (  1 , R  + 1  )  *CGN ST P4-2 1  •  — l./BETA) *PRSXP21  1 , 

1 R  J -PR SXP2 (  1,R-1 )*CCI\iSTM-DELPRM 

_ CALL  ThGMAS _ 

KP  l  =  (Si+  1 

CC  198  1=1, NP1 

ISA  P-R£12.«  I  )  =  SQRT  (  PR  SXP2  K  ? ,  J 
CO  195  1=  1 ,  IMP  1 

IF1 ABS(PRS(2, I )-PRS(3,I ) ) .LE.DELPRS)  GC  TG  197 

-AG  1  GG  194  J=lfMP1 _ _ _ _ _ _ 

PRSXP2I3,J)=PrSXP2(2»J) 

194  PRS(3, J)=PRS(2, J) 

_ _  i-Fl  1  .EQ^NPl  1  GO  TC  113  ______ 

GL  TC  400 

197  IF(I.EQ.NPl)  GC  TC  4C1 

195  CGKTINUH _ 

113  GO  193  J=1,NP1 

PRS(l,J)=PRS<2,J) 

193 PRSXP2I 1, J1=PRS XP2 ( 2 ,  J ) 

KRITE(6,25)  TIP 

25  FORM  AT ( lHO , *  THE  CUMULATIVE  TIME  IS  • , F6 . 3  »  *  M I  MUTES  * ) 

_ laiRlTE(6t27)  (PRSi  1.1  )  ,1  =  1  ,NP1  ) _ 

27  FORMAT ( IX,  14F9.2) 

DC  192  1=1,  NP1 

131 . P&Sfl  ( l  I»PRSJLLU1Z£JLSLI _ _ 

IdiK I  TE  (  6 , 22  )  CPRSQC1  l»I=i,NPl) 

22  FORMAT (1X,14F8.5) 

2CC  CukTIkUE _ _ _ _ 

STOP 

END 
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SGBROU  T  1  ME  SUBEIA 
INTEGfcR  R 

- CJDuMJlOiiZCAL CB/S.U  fc K T  L  (  3  )  ,  PR PG »  n t  TA » B £ T  A  I^uB _ 

COMMON /FUN SUN/ DfcLP, rtGFPRS ,DELPKG,FC HECK 
CGMMGN  PRS(3,105) 

- PKSAV-=l«/fc.*4PRSll.K-1  l+PKSf  1  «  R  )  +  P  R  s  <  1  T  P  »1  )  +  P  >  -S  (  H  T  R-  1  I 

1+PR3(3,K) +  PRS( 3,R  +  1 ) ) 
lF(SLBRTC{2i .LT.O.C)  GG  TG  10 
_ -CALL  FUNPRS 

10  I F { SLLRTC<3) .LT.O.O)  GG  TC  20 
GALL  PRPBLY 

2D  BFTA=BETAI*PRPG»PRSA\/*wnPPKS _ 

RETURN 

END 
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SUBROUTINE  THOMAS 

cCMM0N/MANE/Pk$XP2 { 3 , 105),B(10b,3),D(100),K 
D 1  ME  IVSJjfliSt  hTlQQl tOilQO  j,  6 ( IQ Q 1 
W  (  1  )  =  b (  i  ,2 ) 

CG  10  1=2,  K 

iiL.wm=e(i.2)-B(  i  i-i.jj/At  i-i ) _ 

KM i  =  K—  1 
CC  AO  1=1 » KM  1 
AO  Q  t.Il=8  II  ,3*/M(  I  ) 
u(l)=C(l)/w(  l) 

CC  20  1  =  2, K 

2C  G(I  )  =  LU(  1  )  — d(  1  ,  1)*G(  1-1  ?  )/d  (  1  ) _ 

PKSXP2 { 2  ,  K  +  l ) = G ( K ) 

KN  =  0 


CC  30  1=1, KU 
KM=K  M+ 1 

_ Kf\  =  KN-  1  _ 

30  PRSXP2 (2 , K+l-KM )=G ( K-KM )-G { K-KM ) * PRSXP2 ( 2, K+ 1-KN) 
N  P 1  =  1 3 

_ RETURN  . . _ . , . . . . .  _ 

ENl 
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SCBROUT I  NE  PRMBLY 

COMMON/ PE RM/PECHEK  , DELL TH , DEL PKM  , N 
CGMMGN/JCAL  Cb/bUB  fs.T  C  13)  ,PRMb  ,  BETA,  BET  A1  ,R 
CGMMON/FUNSUN/ DEL P, WOFPRS ,OELPKG, FC HECK 
REAL  LTHC ( 100) , PRMGVT ( 100 ) 

_ INTEGER  Nl.N.k _ 

I  F  (  PECFiEK  .GT  .0  •  1  )  GG  TG  12 
RE AD ( 5 , 2 5  )  NI 

25  FCRMA1  MX.  13) 

WRITE  16*251  M 

R  E AD ( 5 , 2o )  (  LTHQU  )  ,  1=1  ,NI  ) 

26  FCRIYAT(1X.16F7.4) _ 

WRITE(6,26)  (  LTHQ(  1  )  »  1  =  1 ,  N  I  ) 

REAL ( 5,27  )  { PKMGVT (  1 )  ,1  —  1  ,NI  ) 

21  FGRMAT  C .1 X  ,loF b,2) _ _ _ _ 

WR 1TE ( 6 , 2  7 ) (PKMGVT l I  )  ,  1  =  1  ,Ni  ) 

PEChEK=l . 0 

12  PCHECK=FLGAT  (  R—  1  )/Fl  HAKNii _ 

PCFIB  CK  =  P  CH  ECK— 0  *  5/FLG  A T  (N  ) 

DC  18  I = 1 , N I 

- . -  IF ( PChECK .LE .LTHU( 1X1 . GG  LQ  13 

16  CONTINUE 

13  DELK1=PRKGVT  (  1  ) 

LCP=  1 

PChECK=PCHECK+l  •O/FLGAT  {  N  ) 

DG  14  1=1,  M 

_ . . IF  (  PCHECK .XL*  L  I  HQ  (XIX . GG . IQ  15 

14  CONTINUE 

15  U ELK 2= PKMGVT ( I ) 

_ JCP=1 _ _ _ _ 

PRMG=(DELKl+DELK2)/2. 

IF ( FChECK  .LT  .0. 5)  GG  TG  Id 

_ EFMCP.«iLQ«JC£i . GU  _m.  l6 . . . _ _ __ _ 

L EL PRN=( DtLNz-DELKl )/ (2 . *PRHG) *DELPKG 
GC  TC  17 

16  DELPRM=(X _ _ 

17  RETURN 
END 
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SUBRUUTINE  FUNPRS 

CCMMGN/M ANE/  PKSXP2  (  3  »  1 05 )  ,  B (  105 , 3 ) , D ( 100 ) , K 

- £fillM.CtN/M  i  bLSV/pKSV  lb  »TEM,VCHECK,  VISG 

COHHON/FUNSUN/OELPf  HOFPRS  »  L)  ELPKG  t  PCHECK 
CGFMQN/CALCB/SUBRTC (3 ) t PR KG , B t TA , BET A i , R 

_ CuRPCN  P  R  S  (  3  . 1 Q  5 )  _ 

INTEGER  R 
DIMENSION  W{ 31 

.  .13  DC  1C  1  =  1, 3,2 .  _ _ _ 

PRSVIS^PRS (I f K) 

CALL  VISCOS 

_ LC_fa(U  =  l./VISG _ 

W0FPR3=0.  5*1  W(  1  )+v*(  3)  ) 

FChECK-0 • 0 
Cl L  P  =  C  * 

IF< (PRSXP2I ifR-ll+.Oll .GE*PRSXP2( 1VR+1) .AND. ( PRSXP2( 1 , 
1R  +  1  H-.Ol ) .GE.PRSXP2 ( 1  ,R-1 ) )  GC  TG  12 

- DELP=DELP+FkSXP2I  1  ,8+1  )-P8SXP?(  1  . K- 1) _ 

FCHECK= 1 .0 

12  IP  <1 PRSXP2(3,K-1 >  +  .01 ) .GE.PRSXP2 (3,R+1) .AND. { PRSXP2( 3, 

- Ik  +  U  +  .OIJ  .GL.PKSXP21  3,K-m  GC  TO  20 

DELP»0ELP+PRSXP2(3f R+1)-PRSXP2( 3f R-l) 

PCHECK=1 • c 

_ CCPPRS-O. _ _ 

CG  11  I—  1 »  3 , 2 
PRSATM=PR3( 1 fK 1/14.696 

11  UUPPk  b  -  J  C  r  PUS-.O,  PftSATM*889  .  E-Q5*  (  B.958E-Q4+ 

i 1 . 224 E-Ob * 1  PR SA TM~ 1 . ) + i . 1 39 1 E-0 7* ( PRSATM- 1 .  ) * *2 .  )  / ( 1  A  . 
2696J**2. 

_ CELPKG  =  DU  P _ _ _ _ 

DELP=UELP*UCFPRb/ (8.*WGFPRS  i 
GC  TG  21 

20  D  E  L  P  =  0  « _ _ 

21  RETURN 
END 
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SUBKOoT iNt  ViSCCS 

CCMMGN/V ISCSY/PRSV1 S, TEM, VCHECK, V ISG 

- l£feyCH£CK^IlQ>lA  GO  TO  IQ  ,  . 

T  EM- 5 .  /9  *  * ( i  EM-32. ) 

GELTEM=TEM-25. 

- rEMCF  =  DELTGM*4.55£-Cb  _ 

VCHECK- 1.0 

10  PRSATM^PRSV  IS/14.596 

- V  i  Ooi  =  (  i  Q.Q  *  Q  t.LPH.S  AIM— 1  *  U  i  *  (  0 . 0  8  9  5  8  +  1  P  K  S  A  T  M- 1  .  0  )  *  (  Cl .  0  O  0 

10612+0. 399 7E-05*{PRS A TM-i .0 ) ) ) ) *1. 778E-C4 
VISG^V  ISGI  +  TEMCF 
- K  E  f  U  &  N 


END 
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SUBROUTINE  TEkPCL 
CGMMGN/PGLAT £/AN3WER,XR, TCHECK 
DIMENSION  XVI  30  L,  YV  (3  GJ  ,  XLl  3)  ,XL3 )  ,CI1EFF  (5) 
INTEGER  NX  » N Y , NC 
IFITChECK.GT .0.1)  GU  TG  80 
- READ!  5  #  10)  NX.  NY.  Nr.  _ 

10  FORMAT (2X ,314) 

WRITt(6»10)  NX,NY,NC 

RE Al { 5 , i  1  )  i  XV  {  I  )  ,  I~1  ,  NX  ) 

11  FORMAT (IX,11F7.4) 

kR I T  E ( 6 , 11)  (  X  V  (  I  )  ,1  =  1, NX) 

- READ! 5 t  12)  I  Y  V  (  i)«j=l  ,  N  Y  ) _ 

12  FORMAT  (  IX,  11F  7*2  ) 

wR I T  E ( 6 , 12)  I  YV ( I )  ,  1  =  1# NY) 

_  fCH£CK=l.C 

80  DC  30  1  =  2, NX 

1 F  l  XR  .GT • X V ( 1 ) )  GQ  TG  30 
- IFn.Gf.Nf.)  Gf;  TH  7D  _ _ 

IM2=I-2 
DO  29  J=1 , 3 
X  U)  =  X  V(  I  iv  2  +  J 1 

29  Y  (  J  )  =  Y  V  ( IM2+J1 
GO  TG  40 

30  COMINGF _ _ _ 

40  DO  30  1=1,3 

PRCD=  1.0 

DC  49  K=1,J 

IF( I  .EO.K)  GO  TG  49 

PRCD—PRGC^(XR-X(K) )/(X!l)-X(K)i 

AS.  CONTINUE _ _ _ 

50  CGEF F i 1 ) =  PRGD 
ANSU ER  =  0 . 0 
CG  60  1=1,3 

IF L Yi I ) .l 1 . 1 •0E-06)  GC  TC  60 
AN6  W  Ek  =  ANSWER  +CCEFF  (I  )  *  Y  (  I  ) 

60  CONTINUE _ _ _ _ _ 

GO  TC  90 
70  ANSfcER=YV INC) 

9  0  RETURN  __ 
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Figure  25  Transient  Behavior  -  Series  Core  AUS 
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Figure  26  Transient  Behavior  -  Series  Core  ADS 
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Figure  27  Transient  Behavior  -  Series  Core  ADS 
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Figure  28a  Transient  Behavior  -  Berea  Sandstone  2 
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Figure  28b  Transient  Behavior  -  Berea  Sandstone  2 
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Figure  29a  Transient  Behavior  -  Cut  Berea  Sandstone 
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■figure  30  Transient  Behavior  -  Indiana  Limestone 
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Figure  31  Transient  Behavior  -  Indiana  Limestone 
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Figure  32  Relative  Permeability  Ratio  -  Berea  Sandstone  1 
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